スプライン補間1 - 最適化問題としてのスプライン補間の定式化

シリーズ目次

1. スプライン補間1 - 最適化問題としてのスプライン補間の定式化
2. スプライン補間2 - 多次元・高階微分への拡張

1.1 - はじめに

スプライン補間はスプライン曲線とも呼ばれます.
(Wikipedia, スプライン曲線から引用.)

スプライン曲線（スプラインきょくせん、英語: spline curve）とは、スプラインを使用して表現された曲線のこと。スプラインとは区分多項式（区分的に定義された多項式）の事。数学的な背景や曲線あてはめのようなモデルの推定といった側面もあるが、図学や造形デザインで使われることが多い。

本記事の対象は, 関数当てはめとしてのスプライン補間についてです.

スプライン補間は, $(x,y)$に関するデータ点$(s_1,y_1),(s_2,y_2),...$が与えられたときに$y=f(x)$となる$f$を推定する方法の1つです.
スプライン補間一般には, ラグランジュ補間のような高次方程式を使用する方法では問題となる振動を, 区分的に低次多項式を使うことで抑えることができる, という良い性質を持っています.
スプライン補間と書くと通常は, 1次元3次スプライン補間：$f:\mathbb{R}\to\mathbb{R}$を対象にした区分的3次関数による補間を表します.
本記事では,

• 1次元3次スプラインの概要
• 最適化問題としての1次元3次スプライン
• データを必ずしも通らない, 平滑化3次スプライン

について記述します.

1.2 - 1次元3次スプラインの概要

今, $K$個のデータ点$(s_1,y_1),(s_2,y_2),...,(s_K,y_K)$が与えられているとします. ただし, $s_1<\cdots<s_K$と昇順に並べられています. このとき, 閉区間$x:[s_1,s_K]$で定義された3次スプライン$f(x)$は, 次の性質を満たす区分的3次関数です.

• 閉区間$x:[s_1,s_2],[s_2,s_3],...,[s_{K-1},s_K]$において, $f(x)$は3次関数
• $f(x)$は開区間$x:(s_1,s_K)$で$f$の2階微分までが連続
• $y=f(x)$はすべてのデータ点$(s_1,y_1),...,(s_K,y_K)$を通る

このような関数は一意に定まりません. なぜなら, 次のように条件の数が足りないからです.

• $K-1$個の区間があるので, $K-1$個の3次関数で表されている. 3次関数はパラメータ4個で構成されるので, 未知パラメータは$4K-4$個.
• 束縛条件は$4K-6$個. その内訳は
  • 各区間の両端が固定されていることによる$2K-2$個
  • 隣接区間で1階微分が連続であることによる$K-2$個
  • 隣接区間で2階微分が連続であることによる$K-2$個

したがって, 境界条件を定める必要があります. 代表的な境界条件は次のどちらかです.

• $x=s_1,s_K$における$\frac{\mathrm{d}f}{\mathrm{d}x}$を指定(固定境界条件)
• $x=s_1,s_K$において, $\frac{\mathrm{d}^2f}{\mathrm{d}x^2}=0$(自然境界条件)

特に自然境界条件を満たす解を自然3次スプラインと呼びます.

1.3 - 最適化問題の解としての1次元3次スプライン

突然ですが, $x:[s_1,s_K]$を定義域とする$f(x)$に対する, 次の汎関数最適化問題を考えます.

\begin{gather}
\underset{f}{\mathrm{minimize}}\quad \int_{s_1}^{s_K}\mathrm{d}x\,
\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2, \\
\mathrm{subject}\ \mathrm{to}\quad
y_k-f(s_k)=0\quad (k=1,...,K)\quad \mathrm{and}\ \mathrm{BC.}
\end{gather}

ただしBC(Boundary condition, 境界条件)は, $x=s_1,s_K$において$\frac{\mathrm{d}f}{\mathrm{d}x}$を指定する(固定境界条件)かあるいは, $\frac{\mathrm{d}^2f}{\mathrm{d}x^2}=0$(自然境界条件)とします. この最適化問題に対する(適度に滑らかな)解は, $(s_1,y_1),(s_2,y_2),...,(s_K,y_K)$をデータ点とする, 境界条件を満たす1次元3次スプラインとなります.

---

証明のような解説
束縛条件のうち, 積分の境界でない$y_k-f(s_k)=0\quad (k=2,...,K-1)$についてラグランジュの未定乗数法を用いて書き換えると,

E[f](\mu^{K-2})=\int_{s_1}^{s_K}\mathrm{d}x\,\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2-\sum_{k=2}^{K-1}\mu_k(y_k-f(s_k))

に対する$f$と$\mu_2,...,\mu_{K-1}$に関する境界条件付き最適化問題を解けばよいことがわかります. ただし, $\mu^{K-2}$は$(\mu_2,...,\mu_{K-1})$を表すシンボルです. この場合の境界条件は, $\frac{\mathrm{d}f}{\mathrm{d}x}$に関する固定境界条件あるいは$\frac{\mathrm{d}^2f}{\mathrm{d}x^2}=0$に関する自然境界条件と, $y_1-f(s_1)=0$, $y_K-f(s_K)=0$になります.
前式はディラックのデルタ関数を用いて,

E[f](\mu^{K-2})=\int_{s_1}^{s_K}\mathrm{d}x\,\left\{\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2-\sum_{k=2}^{K-1}\mu_k(y_k-f(x))\delta (x-s_k)\right\}

と書き換えられます.
境界条件を満たすような$f(x)$を見つけたとして, そこから微小な擾乱$\varepsilon g(x)$を加えることを考えます. ただし, 任意の$\varepsilon$について$f(x)+\varepsilon g(x)$が境界条件を満たすように$g(x)$を定めます.
このとき,

E[f+\varepsilon g](\mu^{K-2})-E[f](\mu^{K-2})=\varepsilon \int_{s_1}^{s_K}\mathrm{d}x\,
\left\{
2\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\frac{\mathrm{d}^2g}{\mathrm{d}x^2}
+\sum_{k=2}^{K-1}\mu_k g(x)\delta (x-s_k)
\right\}
+O(\varepsilon^2)

と表されます. ここで,

\int_{s_1}^{s_K}\mathrm{d}x\,\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\frac{\mathrm{d}^2g}{\mathrm{d}x^2}
=
\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\frac{\mathrm{d}g}{\mathrm{d}x}\right]_{s_1}^{s_K}
-\left[\frac{\mathrm{d}^3f}{\mathrm{d}x^3}g\right]_{s_1}^{s_K}
+\int_{s_1}^{s_K}\mathrm{d}x\frac{\mathrm{d}^4f}{\mathrm{d}x^4}g

であり, $\frac{\mathrm{d}^3f}{\mathrm{d}x^3}$が有限であるような$f$を考えることとすると, 第1, 2項が0になるので, 結局のところ

E[f+\varepsilon g](\mu^{K-2})-E[f](\mu^{K-2})=\varepsilon \int_{s_1}^{s_K}\mathrm{d}x\,
g(x)\left\{
2\frac{\mathrm{d}^4f}{\mathrm{d}x^4}
+\sum_{k=2}^{K-1}\mu_k \delta (x-s_k)
\right\}
+O(\varepsilon^2)

です. したがって$f(x)$が最適化問題の解となるためには, $x=s_1,...,s_K$を除いて$\frac{\mathrm{d}^4f}{\mathrm{d}x^4}=0$となればよく, $x=s_1,...,s_K$において$\frac{\mathrm{d}^3f}{\mathrm{d}x^3}$がステップ関数程度の特異性を持つことが許されます. そのような解は1次元3次スプラインです.

---

補足
さらに, 自然3次スプラインは次の汎関数最適化問題の解：
$x:(-\infty,+\infty)$を定義域とする$f(x)$に対し,

\begin{gather}
\underset{f}{\mathrm{minimize}}\quad \int_{-\infty}^{+\infty}\mathrm{d}x\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2, \\
\mathrm{subject}\ \mathrm{to}\quad
y_k-f(s_k)=0\quad (k=1,...,K)\quad \mathrm{and}\ \mathrm{BC.}
\end{gather}

区間$x:(-\infty,s_1],[s_K,+\infty)$では$f$は1次関数となります.

---

前述した最適化問題に登場する目的関数:

\int\mathrm{d}x\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2

は, 物理の連続体力学における「梁の曲げのエネルギー」に相当します. すなわち, スプライン補間は

　　　データ点$(s_1,y_1),(s_2,y_2),...$でピン止めされた梁の形状を求める問題

として考えることもできます.

1.4 - データを必ずしも通らない, 平滑化3次スプライン

1次元3次スプラインでは, データ点$(s_1,y_1),...,(s_K,y_K)$を必ず通るような関数$f(x)$を考えてきました. しかし実問題では, 例えば実験データの補間のように, $y$がノイズを含んでいると見なすことが妥当な場合もあると思います. この場合は厳密にデータ点を通る必然性はなく, むしろデータ点を厳密に通っていなくても, より滑らかな関数を採用したいと考えるでしょう. このような場合に適用できるのが, 平滑化3次スプラインです(平滑化スプラインと対比するために, 1次元3次スプラインを単純3次スプラインと呼ぶ場合があります).
平滑化3次スプラインでは, 自然境界条件を考える場合が多く, 本記事でもその設定を採用します.

(自然境界条件を課す)平滑化3次スプラインでは, データ点を通ることを制約条件として与えるのではなく, 二乗誤差$(y_k-f(s_k))^2$を最適化関数に加えます. すなわち,

E[f] = \sum_{k=1}^{K}(y_k-f(s_k))^2 + \lambda 
\int_{-\infty}^{+\infty}\mathrm{d}x\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2

に関する汎関数最適化問題の解が平滑化3次スプラインとなります. $\lambda$は平滑化パラメータと呼ばれ, 正の値をあらかじめ指定しておきます.
$\lambda\to +0$のとき$f(x)$は単純3次スプラインに一致します. また, $\lambda\to +\infty$のとき$f(x)$は1次関数の最小二乗フィッティングに一致します.

平滑化スプラインの計算法

ここでは平滑化スプラインの具体的な計算法には言及しませんが,
https://www.jstage.jst.go.jp/article/jappstat1971/24/1/24_1_27/_pdf/-char/ja
が参考になると思います.

1.5 - おわりに

初めての記事投稿です. ご容赦ください.
平滑化スプラインの具体的な計算法や平滑化パラメータの決定方法について, 今後補強していこうと思います.