自然な 3 次スプラインによる最小二乗近似

最小二乗近似を作成するには、通常、データの近似元となる空間の基底を用意する必要があります。"自然な" 3 次スプラインの空間の例が示すとおり、基底の明示的な作成が常に簡単であるとは限りません。

このセクションでは、明示的な基底が実際には必要ないことを明白にします。なんらかの形式で近似の空間から内挿するための使用可能な手段があれば十分です。このため、Curve Fitting Toolbox™ スプライン関数がベクトル値関数での動作をサポートしているという事実が、非常に重要となります。

このセクションでは、"自然な" 3 次スプラインによる最小二乗近似の次の点について説明します。

問題

特定のブレーク b(1) < ...< b(l+1) をもつ "自然な" 3 次スプラインの空間 S から特定のデータ (x,y) への最小二乗近似を作成するものとします。

一般的な解決

ブレークシーケンス b のすべての "自然な" 3 次スプラインの線形空間 S の基底 (f1,f2,...,fm) がわかれば、c(1)f1+ c(2)f2+ ... + c(m)fm の形式で最小二乗近似が見つかることがわかっています。ここで、ベクトル c は線形システム A*c = y への最小二乗解で、その係数行列は以下によって与えられます。

A(i,j) = fj(x(i)),   i=1:length(x),  j=1:m .

つまり、c = A\y です。

基底マップの必要性

一般的な解では、基底がわかっている必要があるように見えます。ただし、係数シーケンス c を作成するために必要なのは、行列 A を知っていることのみです。このため、"基底マップ"、たとえば、関数 F を用意するだけで十分です。これで、F(c) が特定の重み付け和 c(1)f1+c(2)f2+... +c(m)fm によって与えられるスプラインを返します。そうすることで、j=1:m の場合に、A の j 番目の列を fnval(F(ej),x) として取得できます。ここで、ej は eye(m) (次数 m の単位行列) の j 番目の列です。

さらに好都合なことに、Curve Fitting Toolbox スプライン関数は "ベクトル値" 関数を処理できるため、基底マップ F を作成して、ベクトル値係数 c(i) も同様に処理できます。ただし、取り決めにより、このツールボックス内のベクトル値係数は "列" ベクトルです。したがって、シーケンス c は必然的に列ベクトルの行ベクトル、つまり、"行列" となります。このため、F(eye(m)) は、i 番目の成分が基底要素 fi ( i=1:m) であるベクトル値スプラインです。したがって、データサイトのベクトル x が行ベクトルになると仮定すると、fnval(F(eye(m)),x) は (i,j) エントリが x(j) で i の値となる行列、つまり、求めている行列 A の "転置" となります。一方、指摘したとおり、基底マップ F は係数シーケンス c が行ベクトルになる、つまり、ベクトル A\y の転置になると想定しています。そのため、それに応じてデータ値のベクトル y が行ベクトルになると仮定すると、S からデータ (x,y) への最小二乗近似を次のように取得できます。

F(y/fnval(F(eye(m)),x))

確かに、x と y が任意のベクトル (同じ長さの) になるように準備していた場合は、代わりに以下が使用されます。

F(y(:).'/fnval(F(eye(m)),x(:).'))

"自然な" 3 次スプラインの基底マップ

ブレークシーケンス b(1) < ... < b(l+1) をもつ "自然な" 3 次スプラインの線形空間 S の基底マップ F に必要なことは何でしょうか。この線形空間の次元を m とすると、マップ F は、m ベクトルと S の要素との間に線形の 1 対 1 の対応を設定します。ただし、それはまさしく csape(b, . ,'var') が実行することです。

明確にするため、次の関数 F について考えます。

function s = F(c)
s = csape(b,c,'var');

これは、与えられたベクトル c (b と同じ長さ) の場合、ブレークシーケンス b をもち、b(i) (i=1:l+1) で値 c(i) をとる "一意" の "自然な" 3 次スプラインを提供します。一意性が重要です。ベクトル c と結果のスプライン F(c) との間の対応は、必ず、1 対 1 になります。特に、m は length(b) と等しくなります。その他に、関数 f の点 t での値 f(t) は f に線形に依存します。この一意性によって、F(c) が c に線形に依存することが確保されます (c が fnval(F(c),b) に等しく、可逆線形マップの逆は再び線形マップになるため)。

1 行の解

すべてを組み合わせると、次のコードに辿り着きます。

csape(b,y(:).'/fnval(csape(b,eye(length(b)),'var'),x(:).'),...
'var')

ブレークシーケンス b をもつ "自然な" 3 次スプラインによる最小二乗近似です。

適切な外挿の必要性

コマンドによって MATLAB^® で提供されるデータ、たとえば、国勢調査データで試してみましょう。

load census

さらに、年度 1790:10:1990 を cdate として、値を pop として指定します。ブレークシーケンス 1810:40:1970 を使用します。

b = 1810:40:1970;
s = csape(b, ...
pop(:)'/fnval(csape(b,eye(length(b)),'var'),cdate(:)'),'var');
fnplt(s, [1750,2050],2.2);
hold on
plot(cdate,pop,'or');
hold off

3 つの内部ブレークを使用した "自然な" 3 次スプラインによる最小二乗近似を見ると、与えられた値と共に、結果の近似が太い青線で示されています。

これは、適切な近似のように見えますが、"自然な" 3 次スプラインのようには見えません。再確認すると、"自然な" 3 次スプラインは、最初のブレークの左と最後のブレークの右に対して線形でなければなりません。この近似はいずれの条件も満たしていません。これは、以下の事実に起因します。

与えられたデータへの "自然な" 3 次スプライン内挿は、複数のデータサイトが複数の基本区間にわたる区間をもつ pp 型の csape によって提供されます。一方、基本区間外の pp 型の評価は、pp 型の関連する多項式の最後の区分を使用して、つまり、最大次数の外挿によって、MATLAB ppval または Curve Fitting Toolbox スプライン関数 fnval で実行されます。"自然な" 3 次スプラインの場合は、代わりに 2 次外挿が求められます。つまり、最初のブレークの左に対して、最初のブレークの値と勾配が 3 次スプラインと一致する直線が必要です。このような外挿は、fnxtr によって提供されます。"自然な" 3 次スプラインは最初のブレークでゼロの 2 次導関数をもつため、このような外挿は 3 次となり、3 つの一致条件を満たします。同様に、3 次スプラインの最後のブレークを超えると、最後のブレークの値と勾配がスプラインと一致する直線が必要です。これも、fnxtr によって提供されます。

3 つの内部ブレークを使用した "自然な" 3 次スプラインによる最小二乗近似

正しい 1 行の解

以下の 1 行のコードは、ブレークシーケンス b をもつ "自然な" 3 次スプラインによるデータ (x,y) への正しい最小二乗近似を提供します。

fnxtr(csape(b,y(:).'/ ...
    fnval(fnxtr(csape(b,eye(length(b)),'var')),x(:).'),'var'))

しかし、明らかに、かなり長い行です。

以下のコードは、この正しい公式とプロットを使用しています。3 つの内部ブレークを使用した "自然な" 3 次スプラインによる最小二乗近似に示すように、以前のプロットの上により細い赤線で結果の近似を示しています。

 ss = fnxtr(csape(b,pop(:)'/ ...
  fnval(fnxtr(csape(b,eye(length(b)),'var')),cdate(:)'),'var'));
hold on, fnplt(ss,[1750,2050],1.2,'r'),grid, hold off
legend('incorrect approximation','population', ...
'correct approximation')

3 次スプラインによる最小二乗近似

与えられたブレークシーケンス b をもつすべての 3 次スプラインの空間 S によって近似する場合は、1 行の解が完璧に機能します。このために、Curve Fitting Toolbox スプライン関数を使用する必要すらありません。MATLAB spline を利用できるためです。おわかりのとおり、c は b より 2 つ多いエントリを含むシーケンスで、spline(b,c) はブレークシーケンス b をもつ一意の 3 次スプラインを提供し、これは b(i) (すべての i について) で値 c(i+1) をとり、b(1) で勾配 c(1)、b(end) で勾配 c(end) をとります。したがって、spline(b,.) は S の基底マップです。

その他に、spline(b,c,xi) が、この内挿スプラインの xi で値を提供することがわかっています。最後に、spline がベクトル値データを処理できることがわかっています。したがって、以下の 1 行コードは、ブレークシーケンス b をもつ 3 次スプラインによる、データ (x,y) への最小二乗近似を作成します。

spline(b,y(:)'/spline(b,eye(length(b)),x(:)'))