Eigen线性回归实现：如何无需显式构造全1列完成求解？

核心原理 #

线性回归的常规闭式解需要将特征矩阵拼接全1列求解扩展系数，我们可以将系数拆分为特征权重$w$和偏置项$b$，通过去均值计算避免拼接操作，得到的结果和拼接全1列的方案完全等价：

1. 计算原始特征矩阵的列均值$\bar{X}$与标签向量的均值$\bar{y}$
2. 对特征和标签分别做去均值处理，消去偏置项的影响
3. 仅用去均值后的矩阵求解权重$w$
4. 通过均值回算得到偏置$b$

Eigen代码实现 #

依赖Eigen的懒求值广播机制，全程不会显式构造全1列或全1向量，无额外大内存开销：

#include <Eigen/Dense>
#include <utility>

// 输入参数：
// X: n行d列的原始特征矩阵，无全1列
// y: n维样本标签向量
// 返回值：pair<d维权重向量w, 偏置项b>
std::pair<Eigen::VectorXd, double> linear_regression(const Eigen::MatrixXd& X, const Eigen::VectorXd& y) {
    const int sample_num = X.rows();
    // 列方向求均值，仅生成d维的均值行向量，内存开销极小
    const Eigen::RowVectorXd feature_mean = X.colwise().mean();
    const double label_mean = y.mean();

    // 广播去均值：Eigen内部懒求值，不会显式构造和样本量等长的常量向量
    const Eigen::MatrixXd X_centered = X.rowwise() - feature_mean;
    const Eigen::VectorXd y_centered = y.array() - label_mean;

    // 求解权重，可根据需求替换为QR、SVD等更稳定的分解方式
    const Eigen::VectorXd w = (X_centered.transpose() * X_centered).ldlt().solve(X_centered.transpose() * y_centered);
    // 回算偏置
    const double b = label_mean - feature_mean * w;

    return {w, b};
}

开销说明 #

• 内存开销：仅需要额外存储d维的特征均值向量，不会分配和样本量$n$相关的额外内存，相比拼接全1列的方案节省$n$个浮点数的存储空间，大数据量下优势明显。
• 运行开销：数值计算量和常规拼接方案完全等价，且省去了全1列的内存读写开销，实际运行速度不会低于常规方案，无额外运行时损耗。
• 兼容性：求解部分和常规方案完全一致，可直接替换你原本使用的llt()/ldlt()/householderQr()等分解方法，适配不同的数值稳定性需求。

内容的提问来源于stack exchange，提问作者Svalorzen