はじめに
スプライン曲線を用いるスプライン補完はデータをアップサンプリングする場合や、取得間隔が不均一なデータを均等間隔に補正するような場合に大変有用な数学的手法である。
しかし、スプライン曲線の求解は全データに対して一括して求めることが前提の条件となっている。
全実験データが揃っている場合は問題はないが、実応用の場面ではこれでは少々使い勝手が悪い。
例えば時間進行に伴ってデータが次々と得られるような場合、この時系列データにスプライン補完を適用しようとしても、いつまでもデータが全て揃わないため適用できない。
どこかで区切りをつけスプライン補完を適用することはできるが、それ以降のデータとは繋がりは断たれる。
適用区間を継ぎ足した場合でも移動ウィンドウ法でずらした場合においても、得られるスプライン曲線はそれ以前に得られたものとは異なってしまうからである。
本稿はこのような時間進行に伴ってデータが得られていくようなリアルタイム時系列データに対してもスプライン曲線による補完が実行できる方法を考察し、その適用方法を提案する。
スプライン曲線逐次適用の問題点
スプライン曲線は n 個のデータに対しその間の n-1 区間を各々の多項式で表す。
多項式の次数は 2 次元平面内のデータに対しては3次多項式が用いられる。
これはスプライン曲線を物理的な弾性体と見做した場合、その弾性体曲線内の弾性エネルギーが最も小さくなるような曲線が3次曲線だからである。
その求解は n 個のデータ全部対し一括して求める。
そしてデータの個数が n+1 に変化した時、そのスプライン曲線は n 個のデータの時とは一般に違った曲線になる。
このことから、スプライン曲線を区間を区切って求め繋ぎ合わせたり、移動ウィンドウ方式で適用するデータ点の範囲をずらしながら求める手法を使うことができない。
データ点意外の補完値が範囲を変える毎に変化してしまうからである。
以下でそれを具体例で示す。
次のグラフ（Fig.1）は、4点毎にスプライン曲線を求めてつなげた場合と、10点一括でスプライン曲線を求めたものとを比較したグラフである。
グラフの通り、つなげた場合と一括での場合の曲線は一致していない。
Fig.1 区間毎のスプライン曲線をつなげた場合と全区間一括で求めたスプライン曲線
また、次のグラフ（Fig.2）は４点からなる区間を１点ずつずらしながらスプライン曲線を求めたものと、10点一括でスプライン曲線を求めたものを比較したグラフである。
グラフの曲線が一致していないことがわかる。
Fig.2 4点ずつずらしながら求めたスプライン曲線と一括で求めたスプライン曲線
これらのように、区間を区切って求めたスプイラン曲線は、データ点全部を含めた一括で求めたスプライン曲線と違った曲線となってしまう。
スプライン曲線は「全体の歪みエネルギーが最小になるような曲線である」という物理的な意味があるため、全データ点を用いて一括で求めたスプライン曲線が正しいスプライン曲線だと言える。
しかし、時系列データをリアルタイムで解析する場合など逐次でデータが増えていく場合は、データの取得が停止しない限りデータが全て揃わない。
したがってこのような場合は物理的な意味を考えた正しいスプライン曲線を求めることができないことになる。
つまり、逐次解析にスプライン曲線を用いることができない。
以下では、時事刻々と増えていく時系列データに対して、物理的な意味を考えた場合においても正しいスプライン曲線を逐次で求める方法について考察する。
スプライン曲線の通常求解
ここで通常のスプライン曲線の求解法を確認しておく。
スプライン曲線 S(x) はデータ点の間の多項式関数 S_j(x) を継ぎ足したものである。
このデータ点間の多項式関数 S_j(x) を区分多項式と呼ぶ。
スプライン曲線を求めることは、これら全ての区分多項式の係数を求めることと同値である。
データ点 x_j と x_{j+1} の間にある区間の多項式を S_j(x) とする。
２次元平面においては多項式関数は三次関数を用いる。
この区分多項式は
S_j(x) = a_j + b_j(x - x_j) + c_j(x - x_j)^2 + d_j(x - x_j)^3
で表されるとする。
スプライン曲線で補完する対象のデータ点が n+1 [点] だった場合、間の区間は n [区間] ある。
上記区分多項式の係数は 4 つ あり、区分多項式は区間数と同数なので求める係数の数は 4n [個] になる。
この 4n [個] の係数を求めることで全区間の区分多項式が求まり、スプライン曲線が求まる。
n [区間] を繋ぐデータ点の数は、両端の 2つ を除き、n-2 [点] である。
この n-2 [点] において区分多項式 S_j(x), S_{j+1}(x) が繋がるため、不連続が起きてはならない。
微分値も含めて連続であることを要請するため、以下の条件を定める。
曲線が連続である為の条件
S_{j}(x_{j+1}) = S_{j+1}(x_{j+1}) = y_{j+1}
曲線の微分値（１階、２階）が一致する為の条件
S_{j}'(x_{j+1}) = S_{j+1}'(x_{j+1}) \\
S_{j}''(x_{j+1}) = S_{j+1}''(x_{j+1})
次に両端の２点において曲率を 0 とする。
これは 自然境界条件 と呼ばれる。
両端において曲率を 0 とする境界条件
S_{0}''(x_0) = 0 \\
S_{n-1}''(x_n) = 0
これらの条件によって以下の 4n [個] の連立方程式が得られる。
\scriptsize{
\begin{aligned}
\begin{matrix}
S_0(x_0) & = & a_0 & + & 0 & + & 0 & + & 0 & = & y_0 \\
S_0(x_1) & = & a_0 & + & b_0(x_1 - x_0) & + & c_0(x_1 - x_0)^2 & + & d_0(x_1 - x_0)^3 & = & y_1 \\
S_1(x_1) & = & a_1 & + & 0 & + & 0 & + & 0 & = & y_1 \\
S_1(x_2) & = & a_1 & + & b_1(x_2 - x_1) & + & c_1(x_2 - x_1)^2 & + & d_1(x_2 - x_1)^3 & = & y_2 \\
& \vdots & \\
S_j(x_j) & = & a_j & + & 0 & + & 0 & + & 0 & = & y_j \\
S_j(x_{j+1}) & = & a_j & + & b_j(x_{j+1} - x_j) & + & c_j(x_{j+1} - x_j)^2 & + & d_j(x_{j+1} - x_j)^3 & = & y_{j+1} \\
& \vdots & \\
S_{n-1}(x_{n-1}) & = & a_{n-1} & + & 0 & + & 0 & + & 0 & = & y_{n-1} \\
S_{n-1}(x_{n}) & = & a_{n-1} & + & b_{n-1}(x_{n} - x_{n-1}) & + & c_{n-1}(x_{n} - x_{n-1})^2 & + & d_{n-1}(x_{n} - x_{n-1})^3 & = & y_{n} \\
\end{matrix}
\end{aligned}
}
\scriptsize{
\begin{aligned}
\begin{matrix}
S_0'(x_1) & = & b_0 & + & 2c_0{\Delta x_0} & + & 3d_0{\Delta x_0}^2 & = & b_1 & = & S_1'(x_1) \\
S_1'(x_2) & = & b_1 & + & 2c_1{\Delta x_1} & + & 3d_1{\Delta x_1}^2 & = & b_2 & = & S_2'(x_2) \\
& & & & & \vdots \\
S_j'(x_{j+1}) & = & b_j & + & 2c_j{\Delta x_j} & + & 3d_j{\Delta x_j}^2 & = & b_{j+1} & = & S_{j+1}'(x_{j+1}) \\
& & & & & \vdots \\
S_{n-2}'(x_{n-1}) & = & b_{n-2} & + & 2c_{n-2}{\Delta x_{n-2}} & + & 3d_{n-2}{\Delta x_{n-2}}^2 & = & b_{n-1} & = & S_{n-1}'(x_{n-1}) \\
\end{matrix}
\\
\\
\begin{matrix}
S_0''(x_1) & = & 2c_0 & + & 6d_0{\Delta x_0} & = & 2c_1 & = & S_1''(x_1) \\
S_1''(x_2) & = & 2c_1 & + & 6d_1{\Delta x_1} & = & 2c_2 & = & S_2''(x_2) \\
&&&\vdots \\
S_j''(x_{j+1}) & = & 2c_j & + & 6d_j{\Delta x_j} & = & 2c_{j+1} & = & S_{j+1}''(x_{j+1}) \\
&&&\vdots \\
S_{n-2}''(x_{n-1}) & = & 2c_{n-2} & + & 6d_{n-2}{\Delta x_{n-2}} & = & 2c_{n-1} & = & S_{n-1}''(x_{n-1}) \\
\end{matrix}
\end{aligned}
}
\scriptsize{
\begin{aligned}
\begin{matrix}
S_0''(x_1) & = & 2c_0 & & & = & 0 \\
S_{n-1}''(x_n) & = & 2c_{n-1} & + & 6d_{n-1}{\Delta x_{n-1}} & = & 0 \\
\end{matrix}
\end{aligned}
}
ここで \Delta x_j = (x_{j+1} - x_j) とする。
この連立方程式から以下の行列方程式が作られる。
\tiny{
X=
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & \Delta x_0 & {\Delta x_0}^2 & {\Delta x_0}^3 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & \Delta x_1 & {\Delta x_1}^2 & {\Delta x_1}^3 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 1 & \Delta x_{n-1} & {\Delta x_{n-1}}^2 & {\Delta x_{n-1}}^3 \\
0 & 1 & 2{\Delta x_0} & 3{\Delta x_0}^2 & 0 & -1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 3{\Delta x_0} & 0 & 0 & -1 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 1 & 2{\Delta x_{n-2}} & 3{\Delta x_{n-2}}^2 & 0 & -1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 & 3{\Delta x_{n-2}} & 0 & 0 & -1 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 3{\Delta x_{n-1}} \\
\end{bmatrix}
}
A = \begin{bmatrix}
a_0 \\ b_0 \\ c_0 \\ d_0 \\ a_1 \\ b_1 \\ c_1 \\ d_1 \\ \vdots \\ a_{n-1} \\ b_{n-1} \\ c_{n-1} \\ d_{n-1}
\end{bmatrix}
, \quad
Y = \begin{bmatrix}
y_0 \\ y_1 \\ y_1 \\ y_2 \\ y_2 \\ \vdots \\ y_{n-1} \\ y_{n-1} \\ y_n \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ 0 \\ 0
\end{bmatrix}
この行列Xの逆行列を使い、以下で全ての係数が求まる。
別解
また別解として、条件式を整理すると以下の関係式が導かれる。
\begin{aligned}
a_j &= y_j \\
b_j &= \frac{y_{j+1} - y_j}{\Delta x_j} - \frac{{\Delta x_j}(c_{j+1} + 2c_j)}{3} \\
d_j &= \frac{c_{j+1} - c_j}{3{\Delta x_j}}
\end{aligned}
{\Delta x_j}c_j + 2(\Delta x_{j+1} + \Delta x_j)c_{j+1} + {\Delta x_{j+1}}c_{j+2} = \frac{3(y_{j+2} - y_{j+1})}{\Delta x_{j+1}} - \frac{3(y_{j+1} - y_{j})}{\Delta x_{j}}
式からわかるように、各 c_j を求めれば、残りの a_j, b_j, d_j も導かれる。
各 c_j は
\begin{aligned}
c_0 &= 0 \\
c_{n-1} &= 0
\end{aligned}
であることを用いると、
\tiny{
\begin{aligned}
\begin{bmatrix}
2(\Delta x_1 + \Delta x_0) & \Delta x_1 & & & & & \\
\Delta x_1 & 2(\Delta x_2 + \Delta x_1) & \Delta x_2 & & & & \\
& & & \ddots & \ddots & \ddots & & \text{\huge{0}} & \\
& \text{\huge{0}} & & & \ddots & \ddots & \ddots & & & \\
& & & & & \Delta x_{n-4} & 2(\Delta x_{n-3} + \Delta x_{n-4}) & \Delta x_{n-3} \\
& & & & & & \Delta x_{n-3} & 2(\Delta x_{n-2} + \Delta x_{n-3}) \\
\end{bmatrix}
\begin{bmatrix}
c_1 \\ c_2 \\ \vdots \\ c_{n-2}
\end{bmatrix}
=
\begin{bmatrix}
\frac{3(y_{2} - y_{1})}{\Delta x_{1}} - \frac{3(y_{1} - y_{0})}{\Delta x_{0}} \\
\frac{3(y_{3} - y_{2})}{\Delta x_{2}} - \frac{3(y_{2} - y_{1})}{\Delta x_{1}} \\
\vdots \\
\frac{3(y_{n-2} - y_{n-3})}{\Delta x_{n-3}} - \frac{3(y_{n-3} - y_{n-4})}{\Delta x_{n-4}} \\
\frac{3(y_{n-1} - y_{n-2})}{\Delta x_{n-2}} - \frac{3(y_{n-2} - y_{n-3})}{\Delta x_{n-3}}
\end{bmatrix}
\end{aligned}
}
このような行列方程式になる。
この行列方程式の左辺正方行列は三重対角行列になっているため、トーマスアルゴリズム（Thomas’ algorithm）を用いて解くことができる。
以上がスプライン曲線の通常求解法である。
時間進行による既存区間の曲線変化
データ点が \{p_i\}_{i=0}^{m} 、つまり p_0, p_1, \cdots, p_m の m [点] あるとし、この m [点] のデータ点に対して通常通りにスプライン曲線を求める。
尚、データ点 p_j と p_{j+1} の x 座標、x_j と x_{j+1} は常に x_j < x_{j+1} であるとする。
ここである区間 [p_j, p_{j+1}] に着目し、区間 j とする。
Fig.3 m[点]のデータ点と区間 j
いまデータ点集合 \{p_i\}_{i=0}^{m} に対し新たなデータ点 p_{m+1} が加わったとし、再度 m+1 [点] に対しスプライン曲線を求める。
以後、新たなデータ点が加わるたびにスプライン曲線を更新していくとする。
この時、区間 j の区分曲線はどのような変化をしていくだろうか。
Fig.4 データ点が追加されることで変化するスプライン曲線
区分曲線を決定しているのは区分多項式であるので、区間 j の区分多項式の変化を見る為には 4 つの係数 a_j, b_j, c_j, d_j がどう変化するかに着目すればよい。
Fig.5 は上記の場合に区間 j の係数 a_j, b_j, c_j, d_j がどう変化したかをグラフにしたものの一例である。
Fig.5 区間 j の係数変化の一例
この結果からは、データ点が追加されるにつれ 4 つの係数はいずれも一定値に収束しているようである。
では他のスプライン曲線の場合においても、データ点が追加されるにつれて 4 つの係数はいずれも一定値に収束するのだろうか。
上記を確かめるために数値計算でシミュレーションを実施した。
データ点を一点追加する度に、データ点追加前後の各係数の差を計算してその経過を記録するという手法を用いた。
データ点追加前後の係数値の差が 0 に収束していけば、係数値がある一定値に収束していくことが分かる。
\begin{aligned}
\Delta a_{n} &= a_{j, n+1} - a_{j, n} \\
\Delta b_{n} &= b_{j, n+1} - b_{j, n} \\
\Delta c_{n} &= c_{j, n+1} - c_{j, n} \\
\Delta d_{n} &= d_{j, n+1} - d_{j, n}
\end{aligned}
シミュレーション条件は以下である。
シミュレーション条件
- 区間 j と初期データ点数 m はランダムに設定 (j < m)
- 追加のデータ点 p_{m+n} は、p_{m+1} ～ p_{m+12} の 12 [点]を順に追加（1 \leqq n \leqq 12）
- 追加のデータ点 p_{m+n} の x, y 座標は毎回ランダムに生成
上記、p_{m+1} ～ p_{m+12} の12点を順に追加し各係数の追加前後の差を計算した結果である、
\Delta a_{0} \: \cdots \: \Delta a_{11} \\
\Delta b_{0} \: \cdots \: \Delta b_{11} \\
\Delta c_{0} \: \cdots \: \Delta c_{11} \\
\Delta d_{0} \: \cdots \: \Delta d_{11}
の計 60 [個]の値を求める試行を 1 試行とする。
上記条件の下、計 1,000 [回] 試行した結果をヒートマップにした（Fig.6）。
x 軸はデータ点追加のステップ数（n-1）、y 軸は試行回数（0 ～ 999）、色で示す値の高低はデータ点追加前後の係数値の差である。
Fig.6 データ点追加による係数の変化推移（1000回試行結果）
どの係数もデータ点が 6 点追加されると（x 軸のメモリ 5）、それ以降差が 0 であるように見える。
これは係数の値がある一定値に近づいている事を示す。
実際に試行した結果の値を見てみると、今回のシミュレーションである 12 点のデータ点追加の範囲では、追加前後の係数の差は 0 に収束はしなかった。
従って各係数の値も 12 点追加の範囲では一定値に収束しない。
しかし、データ点の追加が 6 点を超えると 4 つの係数それぞれにおいてその差は 10^{-3} 以下に、7 点を超えると 10^{-4} 以下に収まる結果であった。
なお、定数項である係数 a はその性質上、常に一定値であり収束する（変化しない）。
以上の結果より、データ点が順次追加されていく状況において、最後に追加された点より 7 点以前の区間におけるスプライン曲線はほとんど変化しないということがいえる。(Fig.7)
Fig.7 データ点追加の影響をほぼ受けないスプイラン曲線範囲
以後何度か試行を行ったが、上記範囲を超える結果は得られなかった。
本来は数学的な証明が求められるが、
新規追加点より 7 点以上離れた区間のスプライン曲線は、以後のデータ点の追加による影響をほとんど受けない
という上記結果はおおよそ確からしいと考える。
両端の微分値が既知の区分曲線計算法
前節の結果が正しいとしても、ほとんど変化しないとはいえ 10^{-4} 以下の範囲で各係数値の変化は生じうることには留意しなければならない。
逐次計算によるスプライン曲線の継ぎ足しを行う場合、この係数値の変化を考慮せずに単純継ぎ足しを行うと、これまで同様に継ぎ足し点において不連続が生じてしまう。
Fig.8 単純継ぎ足し時に生じる継ぎ足し点での不連続
では継ぎ足し点で不連続が生じないように継ぎ足しを行うにはどうすればよいだろうか。
ここで、継ぎ足し実行時点においての以下の事実に注目する。
- 継ぎ足し点における1階、2階の微分値は既知である
- 継ぎ足し点の次の点での1階、2階の微分値は、継ぎ足し実行時点での全データ点を対象としたスプライン曲線から求めることができる
もし、継ぎ足し区間両端の微分値から継ぎ足し区間の区分曲線を求めることができれば、継ぎ足しによる不連続が生じないようにすることが可能である。
では、両端の微分値から区分曲線を求める方法を考察しよう。
欠損区間における区分曲線の復元
Fig.9 スプライン曲線の欠損区間
Fig.9 のような場合を考える。
全区間のスプライン曲線が求まっていて、ある区間 j の区分曲線が欠損している場合である。
この区分曲線 S_j を両端の p_j, p_{j+1} の値から復元することを試みる。
まず、p_j, p_{j+1} の y の値は既知である。
次に、区分曲線 S_{j-1}, S_{j+1} が存在することから、点 p_j, p_{j+1} における1階、2階の微分値も既知である。
この計 6 つの値は、この区分曲線を求める際の満たすべき境界条件となる。
区分曲線の多項式（区分多項式） S_j は
S_j(x) = a_j + b_j(x - x_j) + c_j(x - x_j)^2 + d_j(x - x_j)^3
なので、これはすなわち、
\begin{aligned}
&S_j(x_j) &&= y_j \\
&S_{j+1}(x_{j+1}) &&= y_{j+1} \\
&S_j'(x_j) &&= b_{j-1} + 2c_{j-1}(x_j - x_{j-1}) + 3d_{j-1}(x_j - x_{j-1})^2\\
&S_{j+1}'(x_{j+1}) &&= b_{j+1} \\
&S_j''(x_j) &&= 2c_{j-1} + 6d_{j-1}(x_j - x_{j-1})\\
&S_{j+1}''(x_{j+1}) &&= 2c_{j+1} \\
\end{aligned}
を満たす3次多項式を求めることを意味する。
しかし、3次多項式を求める際の未知数は係数 a, b, c, d の 4 つ。
一方、満たすべき境界条件は上記 6 つである。
未知数 4 に対し条件式 6 となっており、これは過剰条件となるため通常の連立方程式では求解できない。
そこで、一旦3次多項式を求めることを放棄しよう。
境界条件が 6 つの場合に、連立方程式で解が一意に求解できる多項式は5次多項式である。
では上記境界条件を満たす5次多項式を求めよう。
この5次多項式を U_j とおく。
U_j(x) = \alpha_j + \beta_j (x - x_j) + \gamma_j (x - x_j)^2 + \delta_j (x - x_j)^3 + \varepsilon_j (x - x_j)^4 + \zeta_j (x - x_j)^5
この U_j に上記境界条件を当てはめると、
\begin{aligned}
\small{
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 \\
1 & \Delta x_{j} & {\Delta x_{j}}^2 & {\Delta x_{j}}^3 & {\Delta x_{j}}^4 & {\Delta x_{j}}^5 \\
0 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & 2{\Delta x_{j}} & 3{\Delta x_{j}}^2 & 4{\Delta x_{j}}^3 & 5{\Delta x_{j}}^4 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 2 & 6{\Delta x_j} & 12{\Delta x_j}^2 & 20{\Delta x_j}^3
\end{bmatrix}
\begin{bmatrix}
\alpha_j \\ \beta_j \\ \gamma_j \\ \delta_j \\ \varepsilon_j \\ \zeta_j
\end{bmatrix}
}
&=
\small{
\begin{bmatrix}
y_j \\
y_{j+1} \\
b_{j-1} + 2c_{j-1}{\Delta x_{j-1}} + 3d_{j-1}{\Delta x_{j-1}}^2 \\
b_{j+1} \\
2c_{j-1} + 6d_{j-1}{\Delta x_{j-1}} \\
2c_{j+1}
\end{bmatrix}
}
\\
XA &= Y
\end{aligned}
という行列方程式になる。
X は特殊な場合を覗いて正則な正方行列であるので逆行列が存在し、
で係数を求めることができる。
こうして求めた5次多項式 U_j(x) は、両端で不連続を起こすことなくなめらかに区間 j-1 と区間 j+1 の間をつなぐ、区間 j の区分曲線となる。
Fig.10 境界条件より復元した5次多項式の区分曲線 U_j(x)
U_j と S_j の関係
境界条件より復元した5次多項式 U_j(x) と、本来ここにあるスプライン曲線の区分多項式である3次多項式の S_j(x) との関係を調べよう。
6 つの境界条件は、本来区間 j-1 と区間 j+1 をなめらかにつないでいた S_j(x) にも当てはまる。
この境界条件を満たさなければ、両端点で不連続となるからである。
従って、この境界条件を満たす 2 つの多項式 S_j(x) と U_j(x) は以下の条件を満たす。
\small{
\begin{align}
&U_j(x_j) &&= S_j(x_j) \qquad \tag{1} \\
&U_j'(x_j) &&= S_j'(x_j) \qquad \tag{2} \\
&U_j''(x_j) &&= S_j''(x_j) \qquad \tag{3} \\
&U_j(x_{j+1}) &&= S_j(x_{j+1}) \qquad \tag{4} \\
&U_j'(x_{j+1}) &&= S_j'(x_{j+1}) \qquad \tag{5} \\
&U_j''(x_{j+1}) &&= S_j''(x_{j+1}) \qquad \tag{6} \\
\end{align}
}
この条件より、2式の間の係数の関係を導出しよう。
まず(1)より、
\small{
\begin{aligned}
\alpha_j + \beta_j (x_j - x_j) + \gamma_j (x_j - x_j)^2 + \delta_j (x_j - x_j)^3 + \varepsilon_j (x_j - x_j)^4 + \zeta_j (x_j - x_j)^5
&=
a_j + b_j (x_j - x_j) + c_j (x_j - x_j)^2 + d_j (x_j - x_j)^3 \\
\alpha_j + \beta_j \cdot 0 + \gamma_j \cdot 0 + \delta_j \cdot 0 + \varepsilon_j \cdot 0 + \zeta_j \cdot 0
&=
a_j + b_j \cdot 0 + c_j \cdot 0 + d_j \cdot 0 \\
\end{aligned}
}
\begin{aligned}
\therefore \alpha_j &= a_j
\end{aligned}
同様に(2)(3)より、
\begin{aligned}
\beta_j &= b_j \\
\gamma_j &= c_j
\end{aligned}
が分かる。
次にこの結果と(4)から、
\small{
\begin{aligned}
\alpha_j + \beta_j{\Delta x_j} + \gamma_j{\Delta x_j}^2 + \delta_j{\Delta x_j}^3 + \varepsilon_j{\Delta x_j}^4 + \zeta_j{\Delta x_j}^5
&=
a_j + b_j{\Delta x_j} + c_j{\Delta x_j}^2 + d_j{\Delta x_j}^3 \\
\delta_j + \varepsilon_j{\Delta x_j} + \zeta_j{\Delta x_j}^2
&=
d_j \\
\\
\therefore \qquad
(\delta_j - d_j) + \epsilon_j{\Delta x_j} + \zeta_j{\Delta x_j}^2
&=
0 \\
\end{aligned}
}
同じく(5)(6)から、
\begin{aligned}
\beta_j + 2\gamma_j{\Delta x_j} + 3\delta_j{\Delta x_j}^2 + 4\varepsilon_j{\Delta x_j}^3 + 5\zeta_j{\Delta x_j}^4
&=
b_j + 2c_j{\Delta x_j} + 3d_j{\Delta x_j}^2 \\
3\delta_j + 4\varepsilon_j{\Delta x_j} + 5\zeta_j{\Delta x_j}^2
&=
3d_j \\
\\
\therefore \qquad
3(\delta_j - d_j) + 4\varepsilon_j{\Delta x_j} + 5\zeta_j{\Delta x_j}^2
&=
0 \\
\\
\\
2\gamma_j + 6\delta_j{\Delta x_j} + 12\varepsilon_j{\Delta x_j}^2 + 20\zeta_j{\Delta x_j}^3
&=
2c_j + 6d_j{\Delta x_j} \\
3\delta_j + 6\varepsilon_j{\Delta x_j} + 10\zeta_j{\Delta x_j}^2
&=
3d_j \\
\\
\therefore \qquad
3(\delta_j - d_j) + 6\varepsilon_j{\Delta x_j} + 10\zeta_j{\Delta x_j}^2
&=
0 \\
\end{aligned}
まとめると、
\begin{bmatrix}
1 & \Delta x_j & {\Delta x_j}^2 \\
3 & 4{\Delta x_j} & 5{\Delta x_j}^2 \\
3 & 6{\Delta x_j} & 10{\Delta x_j}^2
\end{bmatrix}
\begin{bmatrix}
\delta_j - d_j \\
\varepsilon_j \\
\zeta_j
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0 \\ 0
\end{bmatrix}
これを解くと、
\begin{bmatrix}
\delta_j - d_j \\
\varepsilon_j \\
\zeta_j
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0 \\ 0
\end{bmatrix}
以上より、境界条件から以下が導かれる。
\begin{aligned}
\alpha_j &=& a_j \\
\beta_j &=& b_j \\
\gamma_j &=& c_j \\
\delta_j &=& d_j \\
\epsilon_j &=& 0 \\
\zeta_j &=& 0
\end{aligned}
これら2式間の係数の関係は、区間 j 内において変化しない。
2式とも、区間 j 内において係数変化がないからである。
よって、これらの関係式より
\begin{aligned}
U_j(x)
&= \alpha_j + \beta_j(x-x_j) + \gamma_j(x-x_j)^2 + \delta_j(x-x_j)^3 + 0 + 0 \\
&= a_j + b_j(x - x_j) + c_j(x - x_j)^2 + d_j(x - x_j)^3 \\
&= S_j(x) \\
\\
\therefore \quad U_j(x) &= S_j(x)
\end{aligned}
であることが導かれる。
結局、S_j(x) と U_j(x) は同一であることが分かった。
すなわち、区間両端の境界条件を満たす5次多項式を求めると、自動的に 4乗、5乗の項の係数が 0 になり、本来ここに存在していた3次多項式である S_j(x) が導かれることが判明した。
Fig.11 区分曲線S_j(x)と復元曲線U_j(x)の一致
本節の結論は、
両端の微分値が既知である区分曲線 S_j(x) を求めるには、両端の6つの境界条件を満たす5次多項式を求めればよい
というものであった。
逐次スプライン曲線の求解法
ここまでの結果を用いて逐次スプライン曲線を求める方法を構成することができる。
前提として、逐次スプライン曲線は
- データ点が増加しても過去分の曲線は変化しない
- 現時点においてのデータ点すべてを使用して求めたスプライン曲線と、既に求めた過去分の曲線は一致する
ものとする。
まず、
追加されたデータ点より7点以前の区間のスプライン曲線はほとんど変化しない
という結果を信用することにする
この結果は
- 端のデータ点が増えても減っても、7点以上離れていればその影響は曲線にほぼ及ばない。
という事と同義である。
すなわち、
- 16点以上のデータで求めたスプライン曲線の中央区間は、その後ほぼ変化しない
ことを意味する。
移動ウィンドウ法を用いて継ぎ足し区間を求める場合、先頭データ点も移動することに注意しなければならない。
そのため、両側から 7 点以上内側の曲線を使用する必要がある。
Fig.12 データ点追加を行ってもほぼ変化がない中央区間
これを用いて、ほぼ変化がない "中央区間" をデータ点追加毎に逐次で求め継ぎ足して行けばよさそうであるが、先に触れたように 10^{-4} の範囲で係数が変化する以上、単純に継ぎ足すと各データ点で不連続が生じてしまう。
曲線の継ぎ足しを実行する際、始端となるデータ点での y値、１階２階微分値は既に決まっており、違う値にするとこの始端にて不連続が生じる。
一方終端においては未定であり、以後のデータ点の追加において破綻が生じないように値を選択できればよい。
これらを考慮すると、以下のようにすればよいことが分かる。
- 継ぎ足し区間の始端においては、前回の継ぎ足し曲線の終端（今回の始端）の y値、１階２階微分値を用いる
- 継ぎ足し区間の終端においては、次のようにして求めた y値、１階２階微分値を用いる
- 継ぎ足し区間を「16点以上のデータで求めたスプライン曲線の中央区間 」とするようなスプライン曲線を求め、中央区間の終端に該当するデータ点の y値、１階２階微分値を求める
Fig.13 求解済逐次スプライン曲線と16点での中央区間および継足区間
これより始端終端からそれぞれ 3 つの条件、計 6 つが出てくる。
これらを境界条件として前節の方法を用いて区間曲線を求め、継ぎ足し曲線とすればよい。
U_j(x) = \alpha_j + \beta_j (x - x_j) + \gamma_j (x - x_j)^2 + \delta_j (x - x_j)^3 + \varepsilon_j (x - x_j)^4 + \zeta_j (x - x_j)^5
以上が逐次スプライン曲線の求解法である。
以下に実行手順をまとめる。
逐次スプライン曲線求解の実行方法
Step.1
はじめに、先頭から16点以上のデータ点を用いてスプライン曲線を求める。
このデータ点範囲を「ウィンドウ区間」とする。
そして中央区間の終端点（Fig.14 の赤点）での
を取得する。
Fig.14 逐次スプライン求解説明図１
Step.2
次にウィンドウ区間を1点ずらし、Step.1 と同じようにこの新たなウィンドウ区間でスプライン曲線を求め、中央区間の終端点での
を取得する。
Fig.15 では、Step.1 で値を取得した点を青点、今回値を取得した点を赤点で示している。
Fig.15 逐次スプライン求解説明図２
Step.3
Step.1 で取得した3つの値を始端の、Step.2 で取得した3つの値を終端の境界条件として、6つの境界条件の下での5次多項式の係数を求める。
\begin{aligned}
\small{
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 \\
1 & \Delta x_{j} & {\Delta x_{j}}^2 & {\Delta x_{j}}^3 & {\Delta x_{j}}^4 & {\Delta x_{j}}^5 \\
0 & 1 & 0 & 0 & 0 & 0 \\
0 & 1 & 2{\Delta x_{j}} & 3{\Delta x_{j}}^2 & 4{\Delta x_{j}}^3 & 5{\Delta x_{j}}^4 \\
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 2 & 6{\Delta x_j} & 12{\Delta x_j}^2 & 20{\Delta x_j}^3
\end{bmatrix}
\begin{bmatrix}
\alpha \\ \beta \\ \gamma \\ \delta \\ \varepsilon \\ \zeta
\end{bmatrix}
}
&=
\small{
\begin{bmatrix}
y_j \\
y_{j+1} \\
b_{j-1} + 2c_{j-1}{\Delta x_{j-1}} + 3d_{j-1}{\Delta x_{j-1}}^2 \\
b_{j+1} \\
2c_{j-1} + 6d_{j-1}{\Delta x_{j-1}} \\
2c_{j+1}
\end{bmatrix}
}
\\
\\
XA &= Y \\
A &= X^{-1}Y
\end{aligned}
Step.4
求めた係数を使った5次多項式を用いて、継足し区間を構成する（Fig.16 の赤線部分）。
U_j(x) = \alpha + \beta(x-x_j) + \gamma(x-x_j)^2 + \delta(x-x_j)^3 + \varepsilon(x-x_j)^4 + \zeta(x-x_j)^5
尚、求めた係数のうち、4次係数、5次係数は数学的には 0 のはずである。
よって、1次から3次の係数4つだけを用いて3次多項式によって継足し区間を構成してもよい。
U_j(x) = S_j(x) = \alpha + \beta(x-x_j) + \gamma(x-x_j)^2 + \delta(x-x_j)^3
しかし、コンピューターによる数値計算の場合には計算の際の丸め誤差などの影響もあり、\varepsilon, \zeta は 0 でない小さい値になることが多い。
その場合、5次多項式のまま構成する方が精度が良い結果となる。
Fig.16 逐次スプライン求解説明図３
Step.5
再度ウィンドウ区間を1点ずらし、同じように新たなウィンドウ区間でスプライン曲線を求め、中央区間の終端点での
を取得する。
Fig.17 での赤点。青点は Step.2 で値を取得したデータ点である。
Fig.17 逐次スプライン求解説明図４
Step.6
Step.2 で取得した3つの値を始端の、Step.5 で取得した3つの値を終端の境界条件として、新たな6つの境界条件の下での5次多項式の係数を求める。
求めた係数を使い、Step.4 と同じように新しい継ぎ足し区間を構成する。
Fig.18 逐次スプライン求解説明図５
Step. more
以後同様に繰り返し、逐次スプライン曲線を構成する。
Animation.1 逐次スプライン求解実行過程
通常求解と逐次求解の一致精度
ここまでで考察し導出した逐次スプライン曲線の妥当性を検証する。
この逐次スプライン曲線は、本来の全データ点一括で求めた場合のスプライン曲線とどの程度一致するだろうか。
通常のスプライン曲線と、逐次でのスプライン曲線を求めて比較してみよう。
Fig.19 はこの2つの異なる方法で求めたスプライン曲線を同じグラフにプロットしたものである。
Fig.19 逐次スプライン曲線と通常スプライン曲線との比較
ここで、
- 青の実線で示されている逐次スプライン曲線を U(x)
- 赤の破線で示されている通常スプライン曲線を S(x)
とする。
この U(x) と S(x) が理想的には同一範囲において一致していることが望ましい。
しかし実際にはデータ点が追加されるにつれ、追加される前後で通常スプライン曲線 S(x) にはわずかながらに変化が生じる。
問題は、
「ほぼ変化がない」とした「先頭より7点前以前の曲線」の変化がどの程度であるか
である。
次を確認しておこう。
- 逐次で求めた逐次スプライン曲線 U(x) は、データ点が増えても既に求めた分の曲線に変化はない。
- 通常スプライン曲線 S(x) は、データ点が追加されるにつれ曲線全体に変化が生じる。
今回の逐次スプライン曲線求解法の妥当性を検証するためには、ある時点で求められている逐次スプライン曲線と、その同じ範囲の通常スプライン曲線を比較し、一致の度合を検証すればよい。
一致度合の検証方法として、平均二乗誤差（MSE）がある。
これは、2つの曲線が平均どの程度離れているかを数値的に示してくれる。
計算式は以下である。
E = \frac{1}{x_b - x_a} \int_{x_a}^{x_b} (U(x) - S(x))^2 dx
数値計算による評価の場合は、離散版である以下の式を用いる。
E = \frac{1}{n} \sum_{i=0}^{n} \left ( U(x_i) - S(x_i) \right ) ^2
今回はこの MSE を計算し、逐次スプライン曲線の妥当性を検証しよう。
この値 E が小さければ、この逐次スプライン曲線求解法はその精度範囲で妥当と言える。
検証結果
グラフ（Fig.20）は下記条件の下、上式によるMSEを計算しその頻度を示したヒストグラムである。
- 通常スプライン曲線 S(x)、逐次スプライン曲線 U(x) ともに、31点のデータ点区間の曲線
- 2曲線の範囲である31点のデータ点区間を 1000[点]で補完
上記の組を 1000組 用意し、1000個のMSEを計算した。
Fig.20 1000組の逐次、通常スプライン曲線に対するMSEによる比較結果
平均値、中央値は以下であった。
- 平均値：3.732669e-08
- 中央値：3.305561e-08
中央値を採用すれば、
\sqrt{ 3.31 \times 10^{-08} } \fallingdotseq 1.82 \times 10^{-4}
となり、2つのスプライン曲線は 10^{-4} のオーダーで一致する結果となった。
この結果は実使用上においての一致精度として十分である。
これより、
この逐次スプライン曲線を求める方法は、全体一括でスプライン曲線を求める本来の方法と ほとんど同一の曲線を得ることができる
と結論づけてよいであろう。
以上で、本稿で提案した逐次スプライン曲線の求解法の妥当性を示すことができた。
結言
今回示した方法により、逐次でスプライン曲線を求めて補完を行っても、本来のスプライン曲線による補完と同一の結果を得ることができる。
これにより今までの問題点であった、
逐次でデータが追加されていくような時系列データにスプライン補完を適用することができない
状況を打開することができた。
過去から未来分まですべてデータがそろった静的データに対してしか用いることができなかった手法を、リアルタイムな動的データに適用できる可能性を広めることができた。
これはPCやモバイルデバイスなどによるリアルタイム時系列解析などに応用が期待できよう。
何かの一助になれば幸いである。
（※ 間違い等あればご指摘、ご連絡くださいませ）
Discussion