C++の応用 - 最小二乗法

概要

最小二乗法は、データ点と予測モデルとの誤差を最小化する手法として、統計学や機械学習で使用されている手法である。

これは、実測値と予測値の差 (残差) の2乗和を最小にすることである。
2乗を用いる理由は、正と負の誤差を同等に扱えること、および、数学的な扱いが容易になるためである。

線形回帰直線は、次式で表現される。

予測モデル:

$y=ax+b$

yは目的変数、xは説明変数、aは傾き、bは切片である。

残差平方和 (RSS : Residual Sum of Squares) は、次式のように表される。

$RSS=\sum _{i=1}^{n}(y_{i}-(ax_{i}+b))^{2}$

$(y_{i},x_{i})$ は実測データの点である。

この残差平方和を最小化するために、aとbに関する偏微分を0とおく。

${\dfrac {\partial RSS}{\partial a}}=2\sum _{i=1}^{n}(y_{i}-(ax_{i}+b))\cdot (-x_{i})=0$

${\dfrac {\partial RSS}{\partial b}}=2\sum _{i=1}^{n}(y_{i}-(ax_{i}+b))\cdot (-1)=0$

これらの方程式を解くことにより、最適なaとbを求めることができる。
ここで、nはデータ点の数である。

詳細な導出過程は、第4回 - 2変数の回帰分析のページを参照すること。

${\begin{aligned}a&={\dfrac {\sum _{i=1}^{n}{(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}}{\sum _{i=1}^{n}{(x_{i}-{\bar {x}})^{2}}}}\\&={\dfrac {C_{xy}}{\sigma _{x}^{2}}}\\&={\dfrac {\mbox{ x と y の 共 分 散 }}{\mbox{ x の 母 分 散 }}}\end{aligned}}$

${\begin{aligned}b&={\dfrac {\sum _{i=1}^{n}{y_{i}}-a\sum _{i=1}^{n}{x_{i}}}{n}}\\&={\bar {y}}-a{\bar {x}}\end{aligned}}$

最小二乗法は単純な線形回帰だけでなく、多変量解析や非線形回帰等の複雑なモデルにも拡張することが可能である。
例えば、多項式回帰では、予測モデルを $y=ax^{2}+bx+c$ のように設定して、同様の原理で係数を求めることができる。

また、正規方程式を行列形式で表現することにより、効率的に計算することも可能である。
特に大規模なデータセットを扱う場合、数値計算ライブラリを使用して効率的に解を求めることができる。

最小二乗法 (直線近似)

誤差の二乗和を最小化するアルゴリズム

基本的な直線近似では、誤差の二乗和を最小化するアルゴリズムを使用する。
これは、正規方程式を解いて、最適な傾きと切片を求める。

数値精度の向上させるため、Kahan summationによる加算誤差の低減を行うことを推奨する。
また、仕様に応じて許容誤差 (epsilon) を調整することが望ましい。

また、多項式近似や重み付き最小二乗法等も実装することを推奨する。

 #include <vector>
 #include <cmath>
 #include <limits>
 #include <stdexcept>
 
 // 最小二乗法による直線近似を行う関数
 // 戻り値 : 傾きa, 切片b -> pair<double, double>
 std::pair<double, double> leastSquares(const std::vector<double> &x, const std::vector<double> &y)
 {
    // 入力値の検証
    if (x.size() != y.size() || x.empty()) {
       throw std::invalid_argument("Invalid input vectors");
    }
 
    const int n = x.size();
    const double epsilon = std::numeric_limits<double>::epsilon() * 1000;  // 数値誤差の許容範囲
 
    // xの値が全て同じかどうかを確認
    bool all_x_same = true;
    for (int i = 1; i < n; i++) {
       if (std::fabs(x[i] - x[0]) > epsilon) {
          all_x_same = false;
          break;
       }
    }
 
    if (all_x_same) {
       throw std::invalid_argument("All x values are the same - cannot fit a line");
    }
 
    // 合計値の計算
    // Kahan summationアルゴリズムを使用して数値誤差を低減
    double sum_x = 0,
           sum_y = 0;
    double sum_xx = 0,
           sum_xy = 0;
 
    double c_x = 0, c_y = 0;    // 補正項
    double c_xx = 0, c_xy = 0;  // 補正項
 
    // 各種合計値の計算
    for (int i = 0; i < n; i++) {
       // sum_xの計算
       double y_x = x[i] - c_x;
       double t_x = sum_x + y_x;
       c_x = (t_x - sum_x) - y_x;
       sum_x = t_x;
 
       // sum_y の計算
       double y_y = y[i] - c_y;
       double t_y = sum_y + y_y;
       c_y = (t_y - sum_y) - y_y;
       sum_y = t_y;
 
       // sum_xx の計算
       double y_xx = x[i] * x[i] - c_xx;
       double t_xx = sum_xx + y_xx;
       c_xx = (t_xx - sum_xx) - y_xx;
       sum_xx = t_xx;
 
       // sum_xy の計算
       double y_xy = x[i] * y[i] - c_xy;
       double t_xy = sum_xy + y_xy;
       c_xy = (t_xy - sum_xy) - y_xy;
       sum_xy = t_xy;
    }
 
    // 分母の確認
    double denominator = n * sum_xx - sum_x * sum_x;
    if (std::fabs(denominator) < epsilon) {
       throw std::invalid_argument("Denominator too close to zero - numerical instability");
    }
 
    // 傾きaの計算
    // a = (n∑xy - ∑x∑y) / (n∑x² - (∑x)²)
    double a = (n * sum_xy - sum_x * sum_y) / denominator;
 
    // 切片bの計算
    // b = (∑y - a∑x) / n
    double b = (sum_y - a * sum_x) / n;
 
    // 結果の妥当性の確認
    if (!std::isfinite(a) || !std::isfinite(b)) {
       throw std::runtime_error("Calculation resulted in non-finite values");
    }
 
    return {a, b};
 }

 // 使用例
 try {
    // サンプルデータ
    std::vector<double> x = {1.0, 2.0, 3.0, 4.0, 5.0};
    std::vector<double> y = {2.1, 3.8, 6.2, 7.9, 9.3};
 
    auto [a, b] = leastSquares(x, y);  // a: 傾き, b: 切片
                                       // y = ax + b の形で近似直線が得られる
 
    std::cout << "Slope (a) : " << a << std::endl;
    std::cout << "Intercept (b) : " << b << std::endl;
 
    // a, bを使用した処理
    // ...略
 }
 catch (const std::invalid_argument& e) {
    // エラー処理
    std::cerr << "Error: " << e.what() << std::endl;
    return -1;
 }

※注意
入力ベクトルxとyは同じサイズである必要がある。
浮動小数点数の計算であるため、数値誤差に注意が必要となる。

分母がゼロになる可能性 (全てのxが同じ値の場合等) があるため、確認処理を追加する必要がある。

Eigenライブラリ

Eigenライブラリのライセンス

Eigen 3.1.1以降、MPL2ライセンスである。
以前のバージョンでは、全てLGPLv3でライセンスされていた。

現在、いくつかの機能がLGPLでライセンスされたサードパーティのコードに依存していることに注意する。
ただし、constrained_cg等のような機能は、EIGEN_MPL2_ONLYプリプロセッサシンボルを定義してコンパイルすることにより、明示的に無効にすることができる。

Eigenは様々なサードパーティライブラリ (<Eigen/*Support>ヘッダ名で認識可能) のインターフェースクラスを提供している。
これらのライブラリを使用する場合は、そのライブラリのライセンスに注意する必要がある。

Eigenライブラリのインストール (Windows)

vcpkgをインストールする。

git clone https://github.com/Microsoft/vcpkg.git
cd vcpkg

./bootstrap-vcpkg.bat

Eigenライブラリの公式Webサイトにアクセスして、ソースコードをダウンロードする。
ダウンロードしたファイルを解凍する。

Eigenライブラリをインストールする。

cd eigen-<バージョン>
vcpkg install eigen3:x64-windows

Eigenライブラリのインストール (Linux)

• パッケージ管理システムからインストールする場合

# RHEL
sudo dnf install eigen3-devel

# SUSE
sudo zypper install eigen3-devel

• ソースコードからインストールする場合

Eigenライブラリのビルドに必要なライブラリをインストールする。

# RHEL
-

# SUSE
sudo zypper install mpfr-devel adolc-devel fftw3-devel superlu-devel suitesparse-devel

Eigenライブラリの公式Webサイトにアクセスして、ソースコードをダウンロードする。
ダウンロードしたファイルを解凍する。

tar xf eigen-<バージョン>.tar.gz
cd eigen-<バージョン>

Eigenライブラリをビルドおよびインストールする。

mkdir build && cd build

cmake -DCMAKE_BUILD_TYPE=Release \
      -DCMAKE_INSTALL_PREFIX=<Eigenライブラリのインストールディレクトリ> \
      ..

make -j $(nproc)
make install

Eigenライブラリの使用例

行列計算ライブラリであるEigenライブラリを使用することもできる。

Eigenを使用する場合、数値計算の安定性が期待できる。

• QR分解を使用することにより、数値的に安定な解を得られる。
• 行列演算が最適化されており、スループットが良い。

行列演算がEigenライブラリによって抽象化されており、複雑な数値計算がライブラリに委譲されているため、実装の簡潔である。

また、複雑な回帰モデルへの拡張が容易であり、重み付き最小二乗法等への対応も可能である。

 #include <Eigen/Dense>
 #include <vector>
 #include <limits>
 #include <stdexcept>
 
 std::pair<double, double> leastSquares(const std::vector<double> &x, const std::vector<double>& y)
 {
    // 入力値の検証
    if (x.size() != y.size() || x.empty()) {
       throw std::invalid_argument("Invalid input vectors");
    }
 
    const int n = x.size();
    const double epsilon = std::numeric_limits<double>::epsilon() * 1000;
 
    // xの値が全て同じかどうかを確認
    bool all_x_same = true;
    for (int i = 1; i < n; i++) {
       if (std::fabs(x[i] - x[0]) > epsilon) {
          all_x_same = false;
          break;
       }
    }
 
    if (all_x_same) {
       throw std::invalid_argument("All x values are the same - cannot fit a line");
    }
 
    // Eigenのデータ構造に変換
    Eigen::MatrixXd A(n, 2);
    Eigen::VectorXd b(n);
 
    for (int i = 0; i < n; i++) {
       A(i, 0) = x[i];  // 第1列 : x値
       A(i, 1) = 1.0;   // 第2列 : 定数項 (1)
       b(i) = y[i];     // y値
    }
 
    // 最小二乗法による解の計算
    // (A^T * A)^(-1) * A^T * b を計算
    Eigen::Vector2d solution;
    try {
       // QR分解を使用して、数値的に安定な解を取得
       solution = A.colPivHouseholderQr().solve(b);
 
       // 解の品質を確認
       double relative_error = (A * solution - b).norm() / b.norm();
       if (relative_error > epsilon) {
          throw std::runtime_error("Poor fit quality - relative error too large");
       }
    }
    catch (const Eigen::RuntimeError &e) {
       throw std::runtime_error("Failed to compute solution: " + std::string(e.what()));
    }
 
    // 結果の妥当性を確認
    if (!solution.allFinite()) {
       throw std::runtime_error("Calculation resulted in non-finite values");
    }
 
    // 傾き, 切片 {slope, intercept}
    return {solution(0), solution(1)};
 }

 // 使用例
 
 try {
    std::vector<double> x = {1.0, 2.0, 3.0, 4.0, 5.0};
    std::vector<double> y = {2.1, 3.8, 6.2, 7.9, 9.3};
 
    auto [slope, intercept] = leastSquares(x, y);
 
    std::cout << "Fitting results : " << std::endl;
    std::cout << "y = " << slope << "x + " << intercept << std::endl;
 
    // 適合度 (R²) の計算
    Eigen::Map<const Eigen::VectorXd> y_vec(y.data(), y.size());
    double y_mean = y_vec.mean();
    double ss_tot = (y_vec.array() - y_mean).square().sum();
 
     Eigen::VectorXd y_pred(y.size());
     for (size_t i = 0; i < x.size(); i++) {
        y_pred(i) = slope * x[i] + intercept;
     }
 
     double ss_res = (y_vec - y_pred).array().square().sum();
     double r_squared = 1.0 - ss_res / ss_tot;
 
     std::cout << "R² = " << r_squared << std::endl;
 }
 catch (const std::exception &e) {
    std::cerr << "Error : " << e.what() << std::endl;
    return -1;
 }