Javaでのニュートン法による最小二乗法の実装
本記事では、Javaを用いてニュートン法を利用した最小二乗法の実装について解説します。数式はMathJaxを用いて記述し、具体例を多く含めています。以下のリンクから各セクションにジャンプできます。
最小二乗法とは
最小二乗法は、観測データに対してモデルをフィッティングする方法の一つです。具体的には、次のような式で定義されます。
観測データ \( (x_i, y_i) \) に対し、モデル \( f(x; \boldsymbol{\theta}) \) が与えられるとします。このとき、最小二乗法は次の目的関数を最小化します:
\[ S(\boldsymbol{\theta}) = \sum_{i=1}^n \left( y_i – f(x_i; \boldsymbol{\theta}) \right)^2 \]
\( \boldsymbol{\theta} \) はモデルパラメータであり、これを最適化することで誤差を最小化します。
ニュートン法とは
ニュートン法は、非線形方程式の解を反復的に求める方法です。関数 \( g(\boldsymbol{\theta}) \) の根を見つける場合、次の更新式を用います:
\[ \boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k – \mathbf{H}^{-1} \nabla g(\boldsymbol{\theta}_k) \]
ここで、\( \nabla g(\boldsymbol{\theta}) \) は勾配ベクトル、\( \mathbf{H} \) はヘッセ行列(2次微分行列)です。
最小二乗法においては、目的関数の勾配とヘッセ行列を次のように計算します:
- 勾配: \[ \nabla S(\boldsymbol{\theta}) = -2 \sum_{i=1}^n \left( y_i – f(x_i; \boldsymbol{\theta}) \right) \nabla f(x_i; \boldsymbol{\theta}) \]
- ヘッセ行列: \[ \mathbf{H} = 2 \sum_{i=1}^n \nabla f(x_i; \boldsymbol{\theta}) \nabla f(x_i; \boldsymbol{\theta})^\top – 2 \sum_{i=1}^n \left( y_i – f(x_i; \boldsymbol{\theta}) \right) \nabla^2 f(x_i; \boldsymbol{\theta}) \]
Javaによる実装
ここでは、Javaでニュートン法を実装し、最小二乗法を用いてモデルのパラメータを最適化する方法を示します。
// 必要なライブラリのインポート
import java.util.Arrays;
public class NewtonLeastSquares {
// モデル関数 f(x; theta)
public static double model(double x, double[] theta) {
return theta[0] * x + theta[1]; // 例: 線形モデル y = theta0 * x + theta1
}
// 目的関数 S(theta) の勾配を計算
public static double[] gradient(double[][] data, double[] theta) {
double[] grad = new double[theta.length];
Arrays.fill(grad, 0.0);
for (double[] point : data) {
double x = point[0];
double y = point[1];
double error = y - model(x, theta);
grad[0] -= 2 * error * x;
grad[1] -= 2 * error;
}
return grad;
}
// ヘッセ行列を計算
public static double[][] hessian(double[][] data, double[] theta) {
double[][] hess = new double[theta.length][theta.length];
for (double[] point : data) {
double x = point[0];
double y = point[1];
hess[0][0] += 2 * x * x;
hess[0][1] += 2 * x;
hess[1][0] += 2 * x;
hess[1][1] += 2;
}
return hess;
}
// ニュートン法による更新
public static double[] newtonStep(double[][] data, double[] theta) {
double[] grad = gradient(data, theta);
double[][] hess = hessian(data, theta);
double[] step = new double[theta.length];
// ヘッセ行列の逆行列を計算(簡略化のため固定サイズ)
double det = hess[0][0] * hess[1][1] - hess[0][1] * hess[1][0];
double[][] invHess = {
{ hess[1][1] / det, -hess[0][1] / det },
{ -hess[1][0] / det, hess[0][0] / det }
};
for (int i = 0; i < theta.length; i++) {
for (int j = 0; j < theta.length; j++) {
step[i] -= invHess[i][j] * grad[j];
}
}
for (int i = 0; i < theta.length; i++) {
theta[i] += step[i];
}
return theta;
}
public static void main(String[] args) {
// 観測データ (x, y)
double[][] data = {
{1, 2}, {2, 4}, {3, 6}, {4, 8}, {5, 10}
};
// 初期パラメータ
double[] theta = {0.0, 0.0}; // [theta0, theta1]
// 反復計算
for (int i = 0; i < 10; i++) {
theta = newtonStep(data, theta);
System.out.println("Iteration " + (i + 1) + ": " + Arrays.toString(theta));
}
}
}
具体例
上記のコードを実行すると、次のような結果が得られます。これは、線形モデル \( y = \theta_0 x + \theta_1 \) に対してニュートン法を適用した例です。
観測データが \( \{(1, 2), (2, 4), (3, 6), (4, 8), (5, 10)\} \) の場合、最終的なパラメータは次のようになります:
- \( \theta_0 \approx 2.0 \)
- \( \theta_1 \approx 0.0 \)
この結果は、観測データが完全に直線 \( y = 2x \) に一致することを意味します。
まとめ
本記事では、Javaを用いてニュートン法による最小二乗法を実装する方法について解説しました。ニュートン法を利用することで、高速かつ精度の高いパラメータ推定が可能です。ただし、ヘッセ行列の逆行列計算や収束性に注意が必要です。これを応用すれば、より複雑なモデルへのフィッティングも実現できます。