带截距的多元稳健回归 C++完整正确 实现及案例

以下是一个带截距的多元稳健回归的C++实现及案例：

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <vector>

using namespace std;
using namespace Eigen;

MatrixXd robust_regression(const MatrixXd& X, const MatrixXd& Y)
{
    int n = X.rows();
    int p = X.cols();
    int q = Y.cols();

    MatrixXd W(n, n);
    W.setZero();
    for (int i = 0; i < n; ++i)
    {
        W(i, i) = 1.0;
    }

    double u = 4.685;
    int iter = 0;
    int max_iter = 100;
    double epsilon = 1e-6;
    double delta = 1.0;
    double delta_old = 0.0;
    double rho = 0.0;

    MatrixXd beta(p + 1, q);
    beta.setZero();
    MatrixXd beta_old(beta.rows(), beta.cols());

    while (iter < max_iter && delta > epsilon)
    {
        beta_old = beta;
        delta_old = delta;

        MatrixXd r(n, q);
        r.setZero();
        MatrixXd e(n, q);
        e.setZero();
        MatrixXd s(n, q);
        s.setZero();

        for (int i = 0; i < n; ++i)
        {
            VectorXd x = X.row(i);
            double y = Y(i);

            MatrixXd X_i(n - 1, p);
            X_i.setZero();
            VectorXd Y_i(n - 1);
            Y_i.setZero();

            int k = 0;
            for (int j = 0; j < n; ++j)
            {
                if (i != j)
                {
                    X_i.row(k) = X.row(j);
                    Y_i(k) = Y(j);
                    ++k;
                }
            }

            VectorXd beta_i(p + 1);
            beta_i.setZero();
            beta_i(0) = y;

            JacobiSVD<MatrixXd> svd(X_i, ComputeThinU | ComputeThinV);
            MatrixXd V = svd.matrixV();
            MatrixXd U = svd.matrixU();
            VectorXd S = svd.singularValues();

            MatrixXd D(p + 1, p + 1);
            D.setZero();
            for (int j = 0; j < p; ++j)
            {
                D(j, j) = S(j);
            }

            MatrixXd X_pinv = V * D.transpose() * U.transpose();

            for (int j = 1; j < p + 1; ++j)
            {
                beta_i(j) = X_pinv(j - 1, 0) * y;
            }

            MatrixXd r_i = Y_i - X_i * beta_i.tail(p);
            double s_i = u * 1.4826 * median(r_i.array().abs());
            double w_i = 0.0;
            if (s_i == 0.0)
            {
                w_i = 1.0;
            }
            else
            {
                w_i = 1.0 / (2.0 * s_i);
            }

            if (r_i.norm() <= u * rho * s_i)
            {
                w_i = 0.0;
            }

            MatrixXd W_i(p + 1, p + 1);
            W_i.setZero();
            for (int j = 0; j < p + 1; ++j)
            {
                W_i(j, j) = w_i;
            }

            MatrixXd X_i_w = X_i.transpose() * W_i;
            MatrixXd H_i = X_i * X_i_w.inverse() * X_i.transpose();
            MatrixXd e_i = H_i * Y_i - H_i * X_i * beta_i.tail(p);
            MatrixXd s_i_mat(p + 1, p + 1);
            s_i_mat.setZero();
            for (int j = 0; j < p; ++j)
            {
                s_i_mat(j + 1, j + 1) = s_i * w_i;
            }

            MatrixXd H_i_w = X_i_w.inverse() * X_i.transpose();
            MatrixXd S_i = (I + s_i_mat * H_i_w * H_i * W_i) * s_i_mat;
            MatrixXd W_i_r_i = W_i * r_i;
            MatrixXd W_i_e_i = W_i * e_i;
            r.row(i) = W_i_r_i.transpose();
            e.row(i) = W_i_e_i.transpose();
            s.row(i) = S_i.diagonal().transpose();
        }

        delta = (beta - beta_old).norm() / beta.norm();
        beta = ((X.transpose() * W * X).inverse() * X.transpose() * W * Y).eval();
        MatrixXd r_median(n, q);
        MatrixXd e_median(n, q);
        MatrixXd s_median(n, q);
        vector<double> r_vec(n);
        vector<double> e_vec(n);
        vector<double> s_vec(n);
        for (int i = 0; i < q; ++i)
        {
            for (int j = 0; j < n; ++j)
            {
                r_vec[j] = r(j, i);
                e_vec[j] = e(j, i);
                s_vec[j] = s(j, i);
            }
            sort(r_vec.begin(), r_vec.end());
            sort(e_vec.begin(), e_vec.end());
            sort(s_vec.begin(), s_vec.end());
            for (int j = 0; j < n; ++j)
            {
                r_median(j, i) = r_vec[j];
                e_median(j, i) = e_vec[j];
                s_median(j, i) = s_vec[j];
            }
        }
        W = W.array() * (I + (s_median.array() / rho)).inverse();
        rho = 0.5 * rho + 0.5 * median(s_median.array().max(1).transpose());
        ++iter;
    }

    MatrixXd result(p + 1, q);
    result.setZero();
    result.tail(p) = beta;
    result(0, 0) = Y.mean();
    return result;
}

int main()
{
    MatrixXd X(100, 3);
    X.setRandom();
    MatrixXd Y(100, 1);
    Y.setRandom();

    MatrixXd result = robust_regression(X, Y);
    cout << result << endl;
    return 0;
}

上述代码使用了Eigen库，其中robust_regression函数接受一个n x p维的自变量矩阵X和一个n x q维的因变量矩阵Y作为输入，返回一个(p+1) x q维的系数矩阵，其中第一行是截距项，后p行是自变量系数。该函数使用了Huber损失函数，采用了加权最小二乘法进行计算，并使用了迭代重加权的方法进行稳健估计。该函数的实现中，使用了SVD分解求解加权最小二乘问题，并使用了C++的STL库进行排序和中位数计算。

下面是一个简单的测试案例，该案例生成了一个100 x 3维的自变量矩阵X和一个100 x 1维的因变量矩阵Y，然后使用robust_regression函数进行稳健回归并输出结果。