【C++】【Eigen】Gauss-Newton法を実装する
観測データをあるモデルで回帰する場合に、最小二乗法がよく使われます。最小二乗法の解法の一つであるGauss-Newton法について調べて見たので、以下にまとめました。また、C++での実装例も紹介します。
多変数関数のNewton法と最小化への適用
Gauss-Newton法の前に、多変数関数にNewton法を適用する場合を見ていきます。
多変数関数のNewton法
多変数関数
$$
\displaystyle
f(x_{1}, \dots, x_{n}) = 0 \tag{1}
$$
の近似解をNewton法を用いて求めることを考えます。
点$(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n})$における1次近似は、
$$
\displaystyle
g(x_{1}, \dots, x_{n}) = \frac{\partial f}{\partial x_{1}}(x_{1} – x^{(k-1)}_{1}) + \dots + \frac{\partial f}{\partial x_{n}}(x_{n} – x^{(k-1)}_{n}) + f(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n}) \tag{2}
$$
と表せます。今、$g(x^{(k)}_{1}, \dots, x^{(k)}_{n}) = 0$となる点への点$(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n})$からの変化量を$\Delta\boldsymbol{x}^{(k)}=(x^{(k)}_{1} – x^{(k-1)}_{1}, \dots, x^{(k)}_{n} – x^{(k-1)}_{n})$で表すと、(2)は、
$$
\displaystyle
g(x^{(k)}_{1}, \dots, x^{(k)}_{n}) =
\begin{pmatrix}
\frac{\partial f}{\partial x_{1}} && \cdots && \frac{\partial f}{\partial x_{n}}
\end{pmatrix}
\Delta\boldsymbol{x}^{(k)} + f(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n}) = 0 \tag{3}
$$
と表すことができます。
連立方程式への適用
次に、多変数関数の連立方程式(1)
$$
\displaystyle
\left\{
\begin{array}{l}
f_{1}(x_{1}, \dots, x_{n}) = 0 \\
\qquad\quad\vdots \\
f_{m}(x_{1}, \dots, x_{n}) = 0
\end{array}
\right. \tag{4}
$$
の近似解をNewton法を用いて求めることを考えます。
(2)(3)より、(4)の1次近似は
$$
\displaystyle
\begin{pmatrix}
g_{1}(x^{(k)}_{1}, \dots, x^{(k)}_{n}) \\
\vdots \\
g_{m}(x^{(k)}_{1}, \dots, x^{(k)}_{n})
\end{pmatrix} =
\begin{pmatrix}
\frac{\partial f_{1}}{\partial x_{1}} & \cdots & \frac{\partial f_{1}}{\partial x_{n}} \\
\vdots & \ddots & \vdots \\
\frac{\partial f_{m}}{\partial x_{1}} & \cdots & \frac{\partial f_{m}}{\partial x_{n}}
\end{pmatrix} \Delta\boldsymbol{x}^{(k)} +
\begin{pmatrix}
f_{1}(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n}) \\
\vdots \\
f_{m}(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n})
\end{pmatrix} = \boldsymbol{0} \tag{5}
$$
で表せます。ここで、
$$
\displaystyle
\begin{pmatrix}
\frac{\partial f_{1}}{\partial x_{1}} & \cdots & \frac{\partial f_{1}}{\partial x_{n}} \\
\vdots & \ddots & \vdots \\
\frac{\partial f_{m}}{\partial x_{1}} & \cdots & \frac{\partial f_{m}}{\partial x_{n}}
\end{pmatrix} = \boldsymbol{J}
$$
と置くと、(5)の連立方程式は(6)のように変形でき、これを解けば近似解が求まるということになります。
$$
\displaystyle
\boldsymbol{J}\Delta\boldsymbol{x}^{(k)}=-
\begin{pmatrix}
f_{1}(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n}) \\
\vdots \\
f_{m}(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n})
\end{pmatrix} \tag{6}
$$
最小化への適用
ある多変数関数$F(x_{1}, \dots, x_{n})$の最小値（厳密には極値）を求めることを考えます。この場合、各変数で$F$を偏微分し$0$と等しいと置くことで、極値が求められます。
$$
\displaystyle
\left\{
\begin{array}{l}
\frac{\partial F}{\partial x_{1}} = 0 \\
\quad\ \vdots \\
\frac{\partial F}{\partial x_{n}} = 0
\end{array}
\right. \tag{7}
$$
(7)の連立方程式は、(4)と同じ形をしていることがわかります。よって、同様に関数$F$の2階微分を
$$
\displaystyle
\begin{pmatrix}
\frac{\partial^{2} F_{1}}{\partial x^{2}_{1}} & \cdots & \frac{\partial^{2} F_{1}}{\partial x_{1} \partial x_{n}} \\
\vdots & \ddots & \vdots \\
\frac{\partial^{2} F_{m}}{\partial x_{n} \partial x_{1}} & \cdots & \frac{\partial^{2} F_{m}}{\partial x^{2}_{n}}
\end{pmatrix} = \boldsymbol{H}
$$
と置けば、
$$
\displaystyle
\boldsymbol{H}\Delta\boldsymbol{x}^{(k)}=-
\begin{pmatrix}
\frac{\partial F_{1}(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n})}{\partial x_{1}} \\
\vdots \\
\frac{\partial F_{m}(x^{(k-1)}_{1}, \dots, x^{(k-1)}_{n})}{\partial x_{n}}
\end{pmatrix} \tag{8}
$$
これを解くことで、極値の近似解が求まります。
Gauss-Newton法
次の最小二乗問題
$$
S(\boldsymbol{x}, \boldsymbol{a}) = \frac{1}{2}\sum^{n}_{i=1}(f(\boldsymbol{x}_{i}, \boldsymbol{a}) – y_{i})^{2} \tag{9}
$$
について考えます。ここで$\boldsymbol{x} = (x_{1}, \dots, x_{n})$は入力値、$\boldsymbol{y} = (y_{1}, \dots, y_{n})$は観測値、$\boldsymbol{a} = (a_{0}, \dots, a_{m})$は関数のパラメータとします。
$S$を最小にするような、パラメータ$\boldsymbol{a}$を求めるとすると、(7)から
$$
\displaystyle
\left\{
\begin{array}{l}
\frac{\partial S}{\partial a_{0}} = \sum^{n}_{i=1}\Bigl\{\frac{\partial f(\boldsymbol{x}_{i}, \boldsymbol{a})}{\partial a_{0}}(f(\boldsymbol{x}_{i}, \boldsymbol{a}) – y_{i})\Bigr\} = 0 \\
\qquad\qquad\qquad\ \vdots \\
\frac{\partial S}{\partial a_{m}} = \sum^{n}_{i=1}\Bigl\{\frac{\partial f(\boldsymbol{x}_{i}, \boldsymbol{a})}{\partial a_{m}}(f(\boldsymbol{x}_{i}, \boldsymbol{a}) – y_{i})\Bigr\} = 0
\end{array}
\right.
$$
と表されますが、これを行列として表すと(10)のようになります。
$$
\displaystyle
\begin{pmatrix}
\frac{\partial S}{\partial a_{0}} \\
\vdots \\
\frac{\partial S}{\partial a_{m}}
\end{pmatrix} =
\begin{pmatrix}
\frac{\partial f(\boldsymbol{x}_{1}, \boldsymbol{a})}{\partial a_{0}} & \cdots & \frac{\partial f(\boldsymbol{x}_{n}, \boldsymbol{a})}{\partial a_{0}} \\
\vdots & \ddots & \vdots \\
\frac{\partial f(\boldsymbol{x}_{1}, \boldsymbol{a})}{\partial a_{m}} & \cdots & \frac{\partial f(\boldsymbol{x}_{n}, \boldsymbol{a})}{\partial a_{m}}
\end{pmatrix}
\begin{pmatrix}
f(\boldsymbol{x}_{1}, \boldsymbol{a}) – y_{1} \\
\vdots \\
f(\boldsymbol{x}_{n}, \boldsymbol{a}) – y_{n}
\end{pmatrix} \tag{10}
$$
ここで関数$f$のJacobian$\boldsymbol{J}_{f}$, 残差ベクトル$\boldsymbol{r}$
$$
\displaystyle
\boldsymbol{J}_{f} =
\begin{pmatrix}
\frac{\partial f(\boldsymbol{x}_{1}, \boldsymbol{a})}{\partial a_{0}} & \cdots & \frac{\partial f(\boldsymbol{x}_{1}, \boldsymbol{a})}{\partial a_{m}} \\
\vdots & \ddots & \vdots \\
\frac{\partial f(\boldsymbol{x}_{n}, \boldsymbol{a})}{\partial a_{0}} & \cdots & \frac{\partial f(\boldsymbol{x}_{n}, \boldsymbol{a})}{\partial a_{m}}
\end{pmatrix},
\boldsymbol{r} =
\begin{pmatrix}
f(\boldsymbol{x}_{1}, \boldsymbol{a}) – y_{1} \\
\vdots \\
f(\boldsymbol{x}_{n}, \boldsymbol{a}) – y_{n}
\end{pmatrix}
$$
を導入すると、
$$
\displaystyle
\begin{pmatrix}
\frac{\partial S}{\partial a_{0}} \\
\vdots \\
\frac{\partial S}{\partial a_{m}}
\end{pmatrix} = \boldsymbol{J}_{f}^{T} \boldsymbol{r} \tag{11}
$$
となります。次に、$S$の2階微分$\boldsymbol{H}$を考えます。$\boldsymbol{H}$の要素$\frac{\partial^2 S}{\partial a_{k} \partial a_{j}}$について見てみると、
$$
\displaystyle
\begin{split}
\frac{\partial^{2} S}{\partial a_{k}\partial a_{j}}
&= \frac{\partial}{\partial a_{k}}\Bigl\{\sum^{n}_{i=1}\Bigl(\frac{\partial f(\boldsymbol{x}_{i}, \boldsymbol{a})}{\partial a_{j}}(f(\boldsymbol{x}_{i}, \boldsymbol{a}) – y_{i})\Bigr)\Bigr\} \\
&= \sum^{n}_{i=1}\Bigl(\frac{\partial^{2}f(\boldsymbol{x}_{i}, \boldsymbol{a})}{\partial a_{k} \partial a_{j}}(f(\boldsymbol{x}_{i}, \boldsymbol{a}) – y_{i}) + \frac{\partial f(\boldsymbol{x}_{i}, \boldsymbol{a})}{\partial a_{j}}\frac{\partial f(\boldsymbol{x}_{i}, \boldsymbol{a})}{\partial a_{k}}\Bigr)
\end{split} \tag{12}
$$
このようになりますが、関数$f$の2階微分は計算が大変なため、この項を省略したものを$\boldsymbol{H}$の近似として考えることにします。前述のJacobianを用いれば、
$$
\displaystyle
\boldsymbol{H} \approx \boldsymbol{J}_{f}^{T} \boldsymbol{J}_{f} \tag{13}
$$
となります。(11)(13)を(8)に適用するとパラメータ$\boldsymbol{a}$の変化量の近似解が得られます。
$$
\displaystyle
\boldsymbol{J}_{f}^{T} \boldsymbol{J}_{f} \Delta \boldsymbol{a}= -\boldsymbol{J}_{f}^{T} \boldsymbol{r} \tag{14}
$$
初期値からスタートし、(14)を用いて変化量$\Delta \boldsymbol{a}$を求め、パラメータを次々と更新していきます。そして(9)がほぼ変化しなくなったところを極小値とみなすのが、Gauss-Newton法によるアルゴリズムです。
Gauss-Newton法は、原理から最小二乗法にしか適用できないことに注意が必要です。
Eigenでの実装
Eigenを使ったC++の実装例です。
点$(x, y) = (0, -0.9), (1, 1.9), (2, 7.3), (3, 13.8), (4, 23.5)$を2次関数$f(x) = a_{0} + a_{1} x + a_{2} x^{2}$で回帰してみました。
ラムダ式で関数$f$とその偏微分$g$を与え、係数$a$の初期値を適当に与えた後に、GaussNewton()を呼び出して計算します。なお、テンプレートで作成しているので、ラムダ式である必要はなく、operator ()を実装したクラスでも動作します。
GitHubにもソースコードをアップしました。
#include <iostream>
#include <Eigen/Eigen>
//---------------------------------------------
// 関数fのJacobianを計算する
//---------------------------------------------
template <typename Grad>
Eigen::MatrixXd CalcJacobian(
Grad & g, Eigen::VectorXd & a, Eigen::VectorXd & x
)
{
Eigen::MatrixXd ret;
ret.resize(x.size(), a.size());
for (unsigned int i = 0; i < x.size(); i++){
auto grad = g(a, x(i));
for(unsigned int j = 0; j < a.size(); j++){
ret.row(i) = grad;
}
}
return ret;
}
//---------------------------------------------
// Gauss-Newton法の計算
//---------------------------------------------
template <typename Func, typename Grad>
bool GaussNewton(
Func & f, Grad & g, Eigen::VectorXd & a, Eigen::VectorXd & x, Eigen::VectorXd & y,
const double eps, const unsigned int max_loop
)
{
// ループ回数
unsigned int loop = 0;
// 残差二乗和
double s_new, s_old = 0;
while(true){
// 関数fのJacobianの計算
Eigen::MatrixXd jac = CalcJacobian(g, a, x);
// 関数fの計算
Eigen::VectorXd fx(x.size());
for(unsigned int i = 0; i < x.size(); i++){
fx(i) = f(a, x(i));
}
// 残差ベクトルrの計算
Eigen::VectorXd r = fx - y;
// 残差二乗和の計算
s_new = r.squaredNorm();
// 収束判定
if(std::fabs(s_new - s_old) < std::fabs(eps)){
break;
}
s_old = s_new;
//---------------------------
// パラメータaの更新
//---------------------------
// Δaの計算
Eigen::MatrixXd h = jac.transpose() * jac;
Eigen::VectorXd delta_a = h.ldlt().solve(-jac.transpose() * r);
// パラメータaの更新
a += delta_a;
// ループ回数のチェック
if(loop >= max_loop){
break;
}
loop++;
}
return (loop < max_loop);
}
//---------------------------------------------
// main()関数
//---------------------------------------------
int main()
{
Eigen::VectorXd a(3), x(5), y(5);
a << 1, 1, 1;
x << 0, 1, 2, 3, 4;
y << -0.9, 1.9, 7.3, 13.8, 23.5;
// 関数f
auto f = [](Eigen::VectorXd & a, double x){
return a(0) + a(1) * x + a(2) * x * x;
};
// 勾配ベクトルg
auto g = [](Eigen::VectorXd & a, double x){
Eigen::VectorXd ret(3);
ret << 1, x, (x * x);
return ret;
};
if(GaussNewton(f, g, a, x, y, 1.0e-6, 5)){
std::cout << a << std::endl;
}
else{
std::cout << "failure." << std::endl;
}
return 0;
}出力結果
-0.891429
1.81286
1.06429Excelで散布図を作り、2次関数で回帰した結果と比較しました。確かに一致していることがわかります。
まとめ
Gauss-Newton法について説明とC++による実装例を紹介しました。
アルゴリズムについてまとめると以下のようになります。
Gauss-Newton法のアルゴリズム
- パラメータ$a$の初期値を決める
- 残差ベクトル$\boldsymbol{r}$と関数$S$を計算する
- 関数$f$のJacobianを計算する($\boldsymbol{J}_{f})$
- $\boldsymbol{J}_{f}^{T} \boldsymbol{J}_{f} \Delta \boldsymbol{a} =-\boldsymbol{J}_{f}^{T} \boldsymbol{r}$からパラメータの変化量$\Delta \boldsymbol{a}$を求める
- $\boldsymbol{a}$を更新する
- 残差ベクトル$\boldsymbol{r}$と関数$S$を計算する
- 関数$S$の変化量がしきい値以下なら計算終了
- 3に戻る