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ライブラリのインストール (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;
 }