最小二乘法的矩阵推导

最小二乘法的矩阵推导

在一个线性回归模型中，假设：
  $y$ 为一个 n × 1 的目标向量；
  $X$ 为一个n × p的特征变量矩阵；
  $\beta$ 为一个 p × 1 的待估参数向量；
  $\epsilon$ 为一个 n × 1 的误差向量;
则线性回归模型可以表示为：
$$y=X\beta +\epsilon$$
其中，ε为误差向量。
  我们的目标是找到一个向量$\beta$，使得模型的预测值$X\beta$与真实值$y$的差异最小。即最小化误差向量的平方和。误差向量的平方和为：
$$S={\epsilon }^T\epsilon =(y-X\beta {)}^T(y-X\beta )$$
  为了找到使误差最小的β，我们对S进行最小化求解,由于是求最小值，对β求导并令导数为零：
$$\begin{array}{rl}\partial S/\partial \beta & =\frac{\partial (y-X\beta {)}^T(y-X\beta )}{\partial \beta }\\ & =\frac{\partial (y^T-{\beta }^T{X}^T)(y-X\beta )}{\partial \beta }\\ & =\frac{\partial (y^{T}y-y^{T}X\beta -{\beta }^T{X}^{T}y+{\beta }^T{X}^{T}X\beta )}{\partial \beta }\\ & =\frac{\partial (y^{T}y-2{\beta }^T{X}^{T}y+{\beta }^T{X}^{T}X\beta )}{\partial \beta }\\ & =2X^{T}Y+2X^{T}X\beta \\ & =0\end{array}$$
则：$$\beta =(X^{T}X{)}^{-1}{X}^{T}Y$$

c++实现如下：

#include <iostream>
#include <Eigen/Dense> // 使用Eigen库进行矩阵计算

const int M = 5; // 样本数量
const int N = 2; // 变量数量

void leastSquares(const Eigen::MatrixXd& X, const Eigen::VectorXd& y, Eigen::VectorXd& coefficients) {
    Eigen::MatrixXd X_t = X.transpose(); // X'的转置
    Eigen::MatrixXd X_t_X = X_t * X; // X'^T * X'
    Eigen::MatrixXd X_t_X_inv = X_t_X.inverse(); // (X'^T * X')^(-1)
    Eigen::MatrixXd X_t_y = X_t * y; // X'^T * y

    coefficients = X_t_X_inv * X_t_y; // (X'^T * X')^(-1) * X'^T * y
}

int main() {
    // 输入矩阵X
    Eigen::MatrixXd X(M, N);
    X << 1, 2,3, 4, 5, 6,7, 8,9, 10;

    // 输出向量y
    Eigen::VectorXd y(M);
    y << 3,6,8,11,14;

    // 构建设计矩阵X'
    Eigen::MatrixXd X_prime(M, N + 1);
    X_prime << Eigen::MatrixXd::Ones(M, 1), X;

    // 最小二乘法求解
    Eigen::VectorXd coefficients(N + 1);
    leastSquares(X_prime, y, coefficients);
    
    std::cout << "Coefficients: " << std::endl;
    std::cout << coefficients << std::endl;
    return 0;
}