C++ 多元线性回归方程模型 求解（Eigen)

时间：2025-03-06 07:07:11

C++ 多元线性回归方程模型 求解


使用的矩阵运算库是Eigen，相关模块介绍如下：

  1. Core Module：
    Core模块是Eigen最基本和最重要的模块，提供矩阵和向量的基本操作以及一些数学函数。Core模块提供了大量常用的矩阵和向量操作，如矩阵加减乘除、向量点积、范数、转置、共轭等，同时还包含许多数学函数，如反三角函数、指数函数、对数函数等。Core模块内的矩阵和向量操作都可以根据用户需要进行动态大小调整和静态大小设定。

  2. Geometry Module：
    Geometry模块主要提供了用于3D几何计算的各种类和方法，如点、直线、圆、球等。Geometry模块与Core模块紧密集成，通过各种变换和运算将几何图形表示为矩阵或向量，从而提供了一种统一的、高效的表示和计算方式。

  3. LU Module：
    LU模块提供了线性求解器中LU分解的实现，可用于求解线性系统。LU模块主要包含LU分解、列主元LU分解、行主元LU分解等功能。LU模块不仅提供了高效的数值计算，还具有完善的算法稳定性和可靠性，在逆矩阵计算、解线性方程组等领域有广泛应用。

  4. QR Module：
    QR模块提供了QR分解的实现，可用于求解线性系统、奇异值分解等问题。QR模块主要包含Householder变换和Givens变换两种QR分解方法，支持稠密矩阵和稀疏矩阵的QR分解。QR模块在数据压缩、降噪、PCA等领域有着广泛的应用。

  5. SVD Module：
    SVD模块提供了奇异值分解的实现，可用于数据压缩、降噪、PCA等问题。SVD模块提供了JacobiSVD和BDCSVD两种奇异值分解方法，支持稠密矩阵和稀疏矩阵的SVD分解。SVD模块是Eigen中最复杂、最耗时的模块之一，但也是最重要、最基础的模块之一。

  6. Sparse Module：
    Sparse模块提供了稀疏矩阵的相关实现，如存储格式、稀疏矩阵向量乘法等。Sparse模块基于Compressed Sparse Column（CSC）和Compressed Sparse Row（CSR）两种存储格式，提供了高效的稀疏矩阵计算能力。Sparse模块支持各种稀疏矩阵操作，如矩阵加减乘除、向量点积、范数、转置等。

  7. Dense Module：
    Dense模块是Eigen中的一个模块，主要提供了稠密矩阵和向量的相关实现。与Sparse模块相对应，Dense模块中的数据存储方式不考虑矩阵或向量中元素的稀疏性，通过二维数组的形式存储矩阵或一维数组的形式存储向量。Dense模块提供了各种矩阵分解的实现，可用于求解线性系统、奇异值分解等问题。

  8. FFT Module：
    FFT模块提供了快速傅里叶变换（FFT）的实现。FFT模块基于Cooley-Tukey算法，支持一维、二维、三维傅里叶变换，以及实序列和复序列的傅里叶变换。FFT模块在信号处理、频谱分析等领域有着广泛的应用。

  9. Unsupported Module：
    Unsupported模块是Eigen内部项目和实验性质的代码集合，提供了各种尚未成型或正式集成的功能，如GPU加速、Tensor的支持等。Unsupported模块的功能十分丰富，但也具有较高的风险和不稳定性，仅可作为参考。

#include <iostream>
#include <Eigen/Dense>
/**
 * @author @还下着雨ZG
 * @brief 多元线性回归模型的梯度求解函数
 * @param[in] X_train 训练集特征
 * @param[in] y_train 训练集标签lable
 * @param[in] w_in 模型的权重向量
 * @param[in] b_in 模型的偏置
*/
std::pair<Eigen::VectorXd,double> compute_gradient(const Eigen::MatrixXd &X_train,
                                                   const Eigen::VectorXd &y_train,
                                                   const Eigen::VectorXd &w_in,const double &b_in)
{
  int m=X_train.rows();
  int dim_features=X_train.cols(); int dim_w=w_in.size();
  if(dim_features!=dim_w){
    std::cout<<"权重w的维度和训练数据X_train维度不匹配！"<<std::endl;
    return std::make_pair(Eigen::VectorXd::Zero(dim_w),0.0);
  }
  Eigen::VectorXd diff_w(dim_w); double diff_b=0;//存导数值
  for(int i=0;i<m;++i){
    double error=(w_in.dot(X_train.row(i))+b_in-y_train[i]);
    for(int j=0;j<dim_w;++j){
      diff_w[j]=diff_w[j]+error*X_train(i,j);
    }
    diff_b=diff_b+error;
  }
  return std::make_pair(diff_w/m,diff_b/m);
}

/**
 * @author @还下着雨ZG
 * @brief 多元线性回归模型的梯度下降法求解目标函数的最优参数
 * @param[in] X_train 训练集特征
 * @param[in] y_train 训练集标签lable
 * @param[in] learning_rate 学习率
 * @param[in] iterations 迭代次数
 * @param[in] w_initial 初始权重
 * @param[in] b_initial 初始偏置值
*/
std::pair<Eigen::VectorXd,double> gradient_decent_optimization(
                                  const Eigen::MatrixXd &X_train,const Eigen::VectorXd &y_train,
                                  const double &learning_rate,const int iterations,
                                  const Eigen::VectorXd &w_initial,const double &b_initial )
{
  Eigen::VectorXd w(w_initial.size());
  double b=0.0;
  w=w_initial;b=b_initial;
  int m=X_train.rows();
  int n=X_train.cols();
  Eigen::VectorXd diff_w(w_initial.size());
  double diff_b;
  for(int iter=0;iter<iterations;++iter){
    auto diff_wb=compute_gradient(X_train,y_train,w,b);
    diff_w=diff_wb.first;diff_b=diff_wb.second;
    w=w-learning_rate*diff_w;b=b-learning_rate*diff_b;
  }
  return std::make_pair(w,b);
}
int main(){
  using namespace std;
  //定义数据集
  Eigen::Matrix<double,3,4> X_train({{2104., 5, 1, 45}, {1416, 3, 2, 40}, {852, 2, 1, 35}});
  Eigen::Vector3d y_train{460, 232, 178};

  Eigen::Vector4d initial_w =Eigen::Vector4d::Zero();
  double initial_b = 0;
  int iterations = 1000;
  double learning_rate = 5.0e-7;
  //确定优化算法：梯度下降算法
  Eigen::Vector4d w_out;
  double b_out;
  auto wb=gradient_decent_optimization(X_train,y_train,learning_rate,iterations,initial_w,initial_b);
  w_out=wb.first;b_out=wb.second;
  //得到多元线性回归方程的参数w和b
  cout<<"w_out=("<<w_out.transpose()<<")\nb_out="<<b_out<<endl;
  //验证目标函数模型的精度
  for(int i=0;i<y_train.rows();i++){
    cout<<"prediction:"<<w_out.dot(X_train.row(i))+b_out<<"\t"<<"real value:"<<y_train[i]<<endl;
  }


  return 0;
}

相关文章