シリーズ目次

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

2.1 - はじめに

この記事は, スプライン補間に関する第2回です.
前回内容：『スプライン補間1 - 最適化問題としてのスプライン補間の定式化』もあわせてお読みください.

前回は, データ点を必ず通る1次元3次スプラインについての説明をし, それがある種の最適化問題の解であるということについて述べました. さらに, データ点を必ず通るのではなく, 「データ点との二乗誤差」と「曲線の曲げのエネルギー(2階微分の大きさ)」をバランスするような解を採用するという, 平滑化3次元スプラインについて説明しました.

本記事では, それを多次元・高階微分へと拡張したいと思います.

2.2 - 多次元平滑化スプラインに関する最適化問題

$f:\mathbb{R}^d\to \mathbb{R}$, つまり$\mathbf{x}\in\mathbb{R}^d$に対する$y=f(\mathbf{x})$を考えます. そしてこれに対応する$K$個のデータ点$(\mathbf{s}_1,y_1),...,(\mathbf{s}_K,y_K)$が与えられているとします. このとき, $f$を推定する問題について考えてみます.
1次元の曲げのエネルギーは$\left[\frac{\mathrm{d}^2f}{\mathrm{d}x^2}\right]^2$と書かれていましたが, 2次元の場合は

\left[\frac{\partial^2f}{\partial x_1^2}\right]^2
+\left[\frac{\partial^2f}{\partial x_1 \partial x_2}\right]^2
+\left[\frac{\partial^2f}{\partial x_2 \partial x_1}\right]^2
+\left[\frac{\partial^2f}{\partial x_2^2}\right]^2,

$d$次元の場合は

\sum_{\nu_1,\nu_2 = 1,...,d}
\left[\frac{\partial^2f}{\partial x_{\nu_1} \partial x_{\nu_2}}\right]^2

と表されます. さらに微分の階数を2でなく$m$とすると,

\sum_{\nu_1,...,\nu_m = 1,...,d}
\left[\frac{\partial^mf}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}\right]^2

と書かれます($x_1,...,x_d$は$\mathbf{x}$の各座標軸を表しています).
したがって,

E[f]=\sum_{k=1}^{K}(y_k-f(\mathbf{s}_k))^2
+\lambda \int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,
\sum_{\nu_1,...,\nu_m = 1,...,d}
\left[\frac{\partial^mf}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}\right]^2

を目的関数とする最適化問題の解を, 次元が$d$, 階数が$m$における平滑化スプラインとして考えることができます. 積分範囲は$\mathbb{R}^d$全体ですので, 1次元の場合の自然スプラインに対応しています.
これはよく, 薄板スプライン(Thin-plate spline)と呼ばれます.

2.3 - 最適化問題の解法

2.3.1 - 微分方程式

それでは, 前述の最適化問題に関する解法を考えていきましょう. 大前提として, $f(\mathbf{x})$が有限である解にしか興味がないことを心に留めておきましょう.
1次元の場合と同様に, 二乗誤差の部分をディラックのデルタ関数を用いて書き換えます.

\begin{split}
E[f]&=\int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,
\sum_{k=1}^{K}(y_k-f(\mathbf{x}))^2\delta (\mathbf{x}-\mathbf{s}_k)\\
&\quad+\lambda \int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,\sum_{\nu_1,...,\nu_m = 1,...,d}
\left[\frac{\partial^mf}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}\right]^2
\end{split}

$f$に微小な擾乱$\varepsilon g(\mathbf{x})$を加え, $E[f+\varepsilon g]-E[f]$の$\varepsilon$の1次の項を考えると,

\begin{split}
\varepsilon^{-1}(E[f+\varepsilon g]-E[f])&=-2\int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,
\sum_{k=1}^{K}(y_k-f(\mathbf{x}))\delta (\mathbf{x}-\mathbf{s}_k)g(\mathbf{x})\\
&\quad+2\lambda \int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,\sum_{\nu_1,...,\nu_m = 1,...,d}
\frac{\partial^mf}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}
\frac{\partial^mg}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}
\end{split}

第2項について, $m$回部分積分を行うと,

\begin{split}
\varepsilon^{-1}(E[f+\varepsilon g]-E[f])&=-2\int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,
\sum_{k=1}^{K}(y_k-f(\mathbf{x}))\delta (\mathbf{x}-\mathbf{s}_k)g(\mathbf{x})\\
&\quad+2(-1)^m\lambda \int_{\mathbb{R}^d} \mathrm{d}\mathbf{x}\,\sum_{\nu_1,...,\nu_m = 1,...,d}
\frac{\partial^{2m}f}{\partial x_{\nu_1}^2 \cdots \partial x_{\nu_m}^2}g(\mathbf{x})
\end{split}

が得られます. ただし, 目的関数内の$\int_{\mathbb{R}^d} \mathrm{d}\mathbf{x},\left[\frac{\partial^mf}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}\right]^2$が収束するので, 無限遠において$f$の$r$階微分$(r\ge m)$は十分速く0に収束しており, 境界積分は0としました.
したがって,

\sum_{\nu_1,...,\nu_m = 1,...,d}\frac{\partial^{2m}f}{\partial x_{\nu_1}^2 \cdots \partial x_{\nu_m}^2}
=\frac{(-1)^m}{\lambda}\sum_{k=1}^{K}(y_k-f(\mathbf{s}_k))\delta(\mathbf{x}-\mathbf{s}_k)

の形であることが必要です. $d$次元ラプラシアン$\triangle_d$は

\triangle_d=\sum_{\nu = 1}^{d}\frac{\partial^{2}}{\partial x_{\nu}^2}

と書かれるので, $f$が満たすべき微分方程式は

\triangle_d^m f=
\frac{(-1)^m}{\lambda}\sum_{k=1}^{K}(y_k-f(\mathbf{s}_k))\delta(\mathbf{x}-\mathbf{s}_k)

となります.

2.3.2 - スプラインの基底関数

$\triangle_d^m f\propto\delta (\mathbf{x})$となるような$f$の関数形(基底関数)を求めましょう.
まず, 原点からの距離$r$のみに依存する関数形$f(r)$(動径分布関数)を特解として考えます.

\triangle_d f(r)=\frac{\mathrm{d}^2f}{\mathrm{d} r^2}+\frac{d-1}{r}\frac{\mathrm{d}f}{\mathrm{d} r}

ですので(Wikipediaより), $m=1$の場合は

f(r)\propto\begin{cases}
r^{2-d} &(d\ne 2)\\
\log r&(d= 2)
\end{cases}

です.
次に, $r$の正のべき乗に$\triangle_d$を作用させたときの挙動を確認しましょう.

$r^{2k}\rightarrow r^{2k-2}\rightarrow \cdots\rightarrow r^{2}\rightarrow 1 \rightarrow 0$
$r^{2k-1}\rightarrow r^{2k-3}\rightarrow \cdots\rightarrow r\rightarrow r^{-1} \rightarrow r^{-3}\rightarrow \cdots$

と変化します. よって, $d$が奇数のときは$f(r)\propto r^{2m-d}$とすれば, $\triangle_d^m f\propto\delta (\mathbf{x})$を得ることができます.
$d$が偶数の場合は$r$の正のべき乗からは所望の$f(r)$は得られません. そこで次に, $r^{2k}\log r$に$\triangle_d$を作用させたときの挙動を確認しましょう.

r^{2k}\log r\rightarrow r^{2k-2}(\log r+c')\rightarrow \cdots\rightarrow r^{2}(\log r+c'')\rightarrow \log r+c''' \rightarrow r^{-2}\rightarrow \cdots

と変化するので, $d$が偶数のときは$f(r)\propto r^{2m-d}\log r$とすれば, $\triangle_d^m f\propto\delta (\mathbf{x})$を得ることができます.
これらの形式が意味を持つのは, $2m-d\ge 1$の場合に限られます(実は, $2m-d\le 0$の場合は最適化問題に滑らかな解が存在しません). 以降は, $2m-d\ge 1$であるものとして議論を進めます.

最後に, $\triangle_d^m r^{2m-d}$,$\triangle_d^m r^{2m-d}\log r$の係数について議論しましょう.
まず, $m=1$について考えます. $V$を原点中心とする半径$\varepsilon$の$d$次元超級とすると, 一般に

\int_V\mathrm{d}V\,\triangle_d f(r)=\int_{\partial V}\mathrm{d}S\,\mathbf{n}\cdot \mathrm{grad}\,f(r)=S_{d-1} \varepsilon^{d-1}\frac{\mathrm{d} f}{\mathrm{d} r}(\varepsilon)

が成り立ちます. ただし, $S_{d-1}$は$d-1$次元単位球面の表面積です. これを用いて,

\begin{gather}
\triangle_d r^{2-d}=(2-d)S_{d-1}\delta(\mathbf{x}) \\
\triangle_2 \log r=S_{1}\delta(\mathbf{x})
\end{gather}

が得られます. よって, 符号にのみ着目すると

\begin{gather}
\triangle_d r^{2m-d}=(-1)^{\frac{d-1}{2}}g\delta(\mathbf{x})\quad (\text{$d$:奇数}) \\
\triangle_d r^{2m-d}\log r=(-1)^{\frac{d}{2}-1}g\delta(\mathbf{x})\quad (\text{$d$:偶数})
\end{gather}

と書けることがわかります. ただし, $g$は$m,d$に依存する正の定数です.

2.3.3 - スプライン関数の一般解

前述の動径分布関数に関する議論の結果から,

\phi(\mathbf{x})=\begin{cases}
(-1)^{m+\frac{d-1}{2}}r^{2m-d} &(\text{$d$:奇数}) \\
(-1)^{m+\frac{d}{2}-1}r^{2m-d}\log r &(\text{$d$:偶数})
\end{cases}

とすると,

\triangle_d^m \phi(\mathbf{x})=(-1)^m g \delta (\mathbf{x})

となります. この$\phi$を用いて微分方程式の一般解は

f(\mathbf{x})=\sum_{k=1}^K a_k \phi(\mathbf{x}-\mathbf{s}_k)+\text{($x_1,...,x_d$に関する$2m-1$次以下の多項式)}

となります. ただし, もともとの最適化問題において意味のある解になるためには, 任意の$\nu_1,...,\nu_m=1,...,d$について$\frac{\partial^mf}{\partial x_{\nu_1} \cdots \partial x_{\nu_m}}$が無限遠で十分速く0に向かう必要があります.
そのためには, ($x_1,...,x_d$に関する$m-1$次以下の多項式)の各項を$h_1(\mathbf{x}),...,h_C(\mathbf{x})$とすると,

f(\mathbf{x})=\sum_{k=1}^K a_k \phi(\mathbf{x}-\mathbf{s}_k)+\sum_{c=1}^C b_c h_c(\mathbf{x})

かつ, 任意の$c$について

\sum_{k=1}^K a_k h_c(\mathbf{s}_k)=0

が成り立つ必要があります. 後者の条件は$\sum_{k=1}^K a_k \phi(\mathbf{x}-\mathbf{s}_k)$の高次の項が0になるように調整する効果があります.

$h_c(\mathbf{x})$は, 例えば$m=3,d=2$のときは

1,\quad x_1,\quad x_2,\quad x_1^2,\quad x_1x_2,\quad x_2^2

となります.

2.3.4 - 係数の決定

最後に, 微分方程式を満足するように$a_k,b_c$を定めましょう. 考えるべき微分方程式は

\triangle_d^m f=
\frac{(-1)^m}{\lambda}\sum_{k=1}^{K}(y_k-f(\mathbf{s}_k))\delta(\mathbf{x}-\mathbf{s}_k)

です. $f(\mathbf{x})$の表式を代入すると, 任意の$k$について

\lambda g \,a_k=
y_k-\sum_{k'=1}^K a_{k'} \phi(\mathbf{s}_k-\mathbf{s}_{k'})-\sum_{c=1}^C b_c h_c(\mathbf{s}_k)

が成立する必要があることがわかります. さらに整理すると,

\begin{gather}
\sum_{k'=1}^K  (\phi(\mathbf{s}_k-\mathbf{s}_{k'})+\lambda g \,\delta_{kk'})a_{k'}+\sum_{c=1}^C b_c h_c(\mathbf{s}_k)=y_k \quad (k=1,...,K)\\
\sum_{k'=1}^K h_c(\mathbf{s}_{k'})a_{k'}=0 \quad (c=1,...,C)
\end{gather}

を満たす必要があることがわかります. このような$a_1,...,a_K,b_1,...,b_C$は一意に定めることができ, その結果が最適問題の解になります.

2.4 - おわりに

本記事で述べた結果は, 前回の記事で述べた1次元3次平滑化スプラインにも適用することができます. すなわち, 区分的3次式で表す替わりに$|x-s_k|^3$の線形結合で表すことで, 本記事と同一の結果を得ることができます.