RBF 補間 (Radial Basis Function Interpolation) の概要と実装

はじめに

RBF 補間とは

RBF 補間とは、RBF (放射基底関数, radial basis function) を用いて入力となる散布データ (scattered data) を補間することです。または、与えられた散布データを元に非線形な関数をフィッティング (近似) することだと考えることもできます。

どういうときに使えるか

下図のように、N 個の点列 (x_1, y_1),  \ldots ,  (x_N,  y_N) が与えられ、それらの点を通るような滑らかな関数 f を求めたいときに、RBF 補間が使えます。

f:id:yuki-koyama:20140504111710p:plain

また、より高次元 (n 次元とします) の空間における補間も可能です。この場合は N 個の点列 (\mathbf{x}_1,  y_1),  \ldots ,  (\mathbf{x}_N,  y_N) が与えられたときに関数 f: \mathbb{R}^n \rightarrow \mathbb{R} を求めることになります。

どういう特徴があるか

  • RBF 補間はノンパラメトリックなフィッティングの一種です。
  • 与えられた点列を必ず通るような関数を求めることができます。
  • 求まった関数は(多くの場合)滑らかです。
  • ノイズが多く密度の高いデータに対しては過学習 (over-fitting) が起きることがあります。
  • 一度事前計算さえすれば、補間値の評価は高速に行うことができます。

RBF ネットワーク

この記事で紹介する内容は、機械学習データマイニングの分野では RBF ネットワークして紹介されることも多いようです。私の専門分野 (Computer Graphics) では散布データを補間する目的で使われるため RBF 補間と呼ぶことが多いようです。この記事ではこちらの呼び方を使うことにします。

追記 (2015/07/12)

PRML の 6.3 章「RBF ネットワーク」にて以下のように紹介されています(以下引用):

歴史的には,RBFが初めて導入されたのは関数補間 (function interpolation) を正確に行うためであった (Powell, 1987).

引用元:C.M. ビショップ. 元田他監訳. (2008). パターン認識機械学習 (下) ベイズ理論による統計的予測. シュプリンガー・ジャパン株式会社.

RBF 補間という考え方がもともとあり、それをニューラルネットワークとして解釈したものが RBF ネットワークと呼ばれることがあるということだと思います。

RBF 補間の考え方

RBF とは

RBF とは、基本的には「その値が原点からの距離にのみ依存する実数値関数」のことを指し、普通 \phi で表します。

\phi(\mathbf{x}) = \left | \mathbf{x} \right |
\phi(\mathbf{x}) = \exp ( - \beta \left | \mathbf{x} \right | ) (\beta は正の実数)

などは代表的な RBF です。

また、「その値が RBF Center \mathbf{c} からの距離のみに依存する実数値関数」も RBF と呼んだりします。この場合は入力 \mathbf{x} に対して \phi(\mathbf{x},  \mathbf{c}) = \phi(\left | \mathbf{x} - \mathbf{c} \right |) となります。

関数の表現

RBF 補間では、ある重み \mathbf{w} = (w_1,  \ldots ,  w_N) を用いて、関数 f

f(\mathbf{x}) = \sum_i^N w_i \phi(\mathbf{x},  \mathbf{x_i})

と表します。つまり、RBF の重み付き線形和です。あとは与えられたデータに対してこの重みを予め上手く見つけておけば良いということになります。

重みの決め方

j 番目の入力データ (\mathbf{x}_j,  y_j) に関して、重み \mathbf{w} が満たすべき条件式として次のような式を考えます。

y_j = \sum_i^N w_i \phi(\mathbf{x}_j,  \mathbf{x}_i)

つまり、得られる関数 f が点 (\mathbf{x}_j,  y_j) を通るという条件です。これを j=1, \ldots , N の全てについて考えると、以下のように列挙できます。

y_1 = \sum_i^N w_i \phi_{1, i}
\vdots
y_N = \sum_i^N w_i \phi_{N, i}

ここで簡単のため \phi_{j, i} = \phi(\mathbf{x}_j,  \mathbf{x}_i) という表記を用いました。すると、これらは行列形式でまとめて

\mathbf{y} = \Phi \mathbf{w}

と書くことができます。ただし \mathbf{y}\Phi はそれぞれ

\mathbf{y} = (y_1,  \ldots ,  y_N),
\Phi = \begin{pmatrix} 
\phi_{1, 1} & \cdots & \phi_{1, N} \\
\vdots & \ddots & \vdots \\
\phi_{N, 1} & \cdots & \phi_{N, N}
\end{pmatrix}

を表しているとします。ちなみに \phi_{j, i} = \phi_{i, j} ですから、\Phi は対称行列になっています。さて、この条件式は \mathbf{w} に関する連立一次方程式と見做すことができ、またほとんどの場合は \Phi には逆行列が存在するので、

\mathbf{w} = \Phi^{-1} \mathbf{y}

とすることで解を求めることができます。

C++, Eigen を用いた RBF 補間の実装

C++ による実装例です。行列演算の部分は、ここでは外部ライブラリとして Eigen を用いることにします。

Eigen

重みの計算方法

入力 (\mathbf{x}_1,  y_1),  \ldots ,  (\mathbf{x}_N,  y_N) に対して重み \mathbf{w} を計算する関数は例えば以下のように書けます。

#include <Eigen/Core>
using namespace Eigen;
using namespace std;
VectorXd computeWeights(VectorXd y, vector<VectorXd> x)
{
    MatrixXd Phi = MatrixXd::Zero(N, N);
    for (int i = 0; i < N; ++ i) {
        for (int j = 0; j < N; ++ j) {
            Phi(i, j) = phi(x[i], x[j]);
        }
    }
    VectorXd w = solveLinearSystem(Phi, y);
    return w;
}

ここで、RBF を表す関数 phi を以下のように定義しておきます。ここでは \phi(\mathbf{x})=\left|\mathbf{x}\right| の定義を用いることにすると、

double phi(VectorXd x)
{
    return x.norm();
}
double phi(VectorXd x, VectorXd c)
{
    return phi(x - c);
}

のようになります。

また一般的な連立一次方程式 \mathbf{A}\mathbf{x}=\mathbf{b} を解くための関数 solveLinearSystem は、どのアルゴリズムを用いても良いですが、例えば完全ピボット選択の LU 分解 (LU decomposition with full pivoting) を用いる場合は

#include <Eigen/LU>
VectorXd solveLinearSystem(MatrixXd A, VectorXd b)
{
    FullPivLU<MatrixXd> lu(A);
    VectorXd x = lu.solve(b);
    return x;
}

のように書きます。今回 \Phi は正定値エルミート行列ですので、コレスキー分解 (Cholesky decomposition) を用いても良いです。一回だけ計算すれば良いとはいえ、この部分の計算量 (時間計算量も空間計算量も) は大きく、例えば N=10000 などになるとノート PC で計算するのが困難になる場合もあるので注意が必要です。

補間結果の値の評価

以上のように求めた重み \mathbf{w} を用いると、任意の入力 \mathbf{x}_{\mathrm{input}} に対する補間結果  y_{\mathrm{interp}} = f(\mathbf{x}_{\mathrm{input}}) は以下のように計算できます。

double f(VectorXd xInput)
{
    double yInterp = 0.0;
    for (int i = 0; i < N; ++ i) {
        yInterp += w[i] * phi(xInput, x[i]);
    }
    return yInterp;
}

ここからも分かる通り、補間値の評価の処理は RBF Center の数 N に対して計算量は O(N) となっています。一度重みの計算さえしておけば、補間値の評価は N が多少大きくても非常に高速に計算可能です。

ダウンロード

以下に C++, Eigen を用いて RBF 補間の機能を実装した例をレポジトリごと公開しました。

yukikoyama / RBF Interpolator — Bitbucket

コンパイルには Qt が必要です。GUI 上で二次元の空間に RBF Center を対話的に追加し、補間の様子をヒートマップで描画するようなプログラムになっています。

上記の本文中で紹介したソースコードとは名前の付け方や型が微妙に違うので注意して下さい。

RBF 補間の応用

Computer Graphics (CG) 分野での応用例

以下の 2001 年の SIGGRAPH で発表された論文では

  1. 物体表面を表す三次元点群データから滑らかな表面の構築
  2. 一部が欠損している表面メッシュデータの穴埋め

の二つの目的を、RBF 補間を応用することで達成しています。

この記事を書いている現在において被引用数が 1300 を超える有名論文で、RBF 補間について詳しく議論しています。Computer Graphics 分野で RBF 補間を扱うときはこの論文を引用するのが便利だと思います。

J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum, and T. R. Evans.
Reconstruction and Representation of 3D Objects with Radial Basis Functions
SIGGRAPH '01
Reconstruction and representation of 3D objects with radial basis functions

また Bickel らによる以下の二本の論文では、次のような処理が行われています。

  1. Measurement-Based Model Fitting. 特定の条件下での実験観測データから、その条件下での物理特性を推定する
  2. Data-Driven Simulation. 推定された特定の条件下における物理特性を RBF Center と見做し、RBF 補間することで、未知の条件下での物理特性を推定し、そのような未知の条件下でも信頼度の高いシミュレーションを行う

少々分かりにくいですが、後半で scattered data interpolation の考え方が使われており、ここでは RBF 補間が採用されています。

B. Bickel, M. Bächer, M. A. Otaduy, W. Matusik, H. Pfister, and M. Gross.
Capture and Modeling of Non-Linear Heterogeneous Soft Tissue
SIGGRAPH '09

B. Bickel, M. Bächer, M. A. Otaduy, H. R. Lee, H. Pfister, M. Gross, and W. Matusik.
Design and Fabrication of Materials with Desired Deformation Behavior
SIGGRAPH '10