Eigen基础用法

基础

定义一个 3∗3 的浮点型矩阵

typedef Eigen::Matrix<float, 3, 3> MyMatrix33f;

定义一个 3∗1 的浮点型列向量

typedef Eigen::Matrix<float, 3, 1> MyVector3f;

定义一个具有动态维数的矩阵

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MyMatrix;

初始化对象

MyMatrix33f a; 
MyVector3f v; 
MyMatrix m(10,15);

初始化矩阵为0矩阵

a = MyMatrix33f::Zero();
std::cout << "Zero matrix:\n" << a << std::endl;

初始化矩阵为单位矩阵

a = MyMatrix33f::Identity();
std::cout << "Identity matrix:\n" << a << std::endl;

添加图片注释，不超过 140 字（可选）

初始化向量中的元素为随机数

v = MyVector3f::Random();
std::cout << "Random vector:\n" << v << std::endl;

添加图片注释，不超过 140 字（可选）

赋值初始化矩阵

v = MyVector3f::Random();
std::cout << "Random vector:\n" << v << std::endl;

添加图片注释，不超过 140 字（可选）

矩阵元素赋值

a(0, 0) = 3;
std::cout << "Matrix with changed element[0][0]:\n" << a << std::endl;

添加图片注释，不超过 140 字（可选）

可以使用Map类型的对象来封装Matrix类型对象中存在的C++数组或向量。这种映射对象会使用来自底层对象的内存和值，而不会分配额外的内存和拷贝值。

int data[] = {1, 2, 3, 4};
Eigen::Map<Eigen::RowVectorXi> v_map(data, 4);
std::cout << "Row vector mapped to array:\n" << v_map << std::endl;
std::cout << "where is the data in memory:\n"<<&data[0] <<std::endl;
std::cout << "where is the v_map in memory:\n"<<v_map.data() <<std::endl;

添加图片注释，不超过 140 字（可选）

矩阵也是一样的

std::vector<float> vdata = {1, 2, 3, 4, 5, 6, 7, 8, 9};
Eigen::Map<MyMatrix33f> a_map(vdata.data());
std::cout << "Matrix mapped to array:\n" << a_map << std::endl;
std::cout << "where is the data in memory:\n"<<&vdata[0] <<std::endl;
std::cout << "where is the a_map in memory:\n"<<a_map.data() <<std::endl;

添加图片注释，不超过 140 字（可选）

Eigen库中的矩阵和向量算术运算可以通过重载标准的C++算术运算符如+、-或*，或者通过方法如dot()和cross()提供。下面的代码示例说明了如何在Eigen中表达一般的数学运算：

定义 2∗2 矩阵

Eigen::Matrix2d a;
a << 1, 2, 3, 4;
Eigen::Matrix2d b;
b << 1, 2, 3, 4;

逐元素相乘/相除/乘4

Eigen::Matrix2d result = a.array() * b.array();
std::cout << "element wise a * b :\n" << result << std::endl;

result = a.array() / b.array();
std::cout << "element wise a / b :\n" << result << std::endl;

a = b.array() * 4;
std::cout << "element wise a = b * 4 :\n" << a << std::endl;

添加图片注释，不超过 140 字（可选）

矩阵加法/乘法

result = a + b;
std::cout << "matrices a + b :\n" << result << std::endl;

a += b;
std::cout << "matrices a += b :\na=\n" << a << std::endl;
std::cout <<"b=\n"<<b<<std::endl;

result = a * b;
std::cout << "matrices a * b :\n" << result << std::endl;

添加图片注释，不超过 140 字（可选）

Eigen的延迟计算机制

在 Eigen 中，像 + 这样的算术运算符不会立刻执行计算，而是返回一个表达式对象，表示“应该怎么计算”，真正的计算是在赋值（=）等操作发生时才进行。

auto expr = a + b;

这时并没有真正把 a 和 b 相加。expr 是一个“表达式对象”，记录的是“加法这个动作”，并不包含结果矩阵。

所以 expr 其实是一个类似“配方”的对象，它告诉 Eigen：“我想要把 a 和 b 相加”

MatrixXd result = expr;

这个时候，Eigen 才会真正执行 a + b 的操作，并把结果写入 result。

但如果你这么写,就不会有延迟计算的问题。

MatrixXd temp = a + b;  // 立刻求值并存储结果

部分运算

随机生成一个 4∗4 的矩阵

Eigen::MatrixXf m = Eigen::MatrixXf::Random(4, 4);
std::cout << "Random 4x4 matrix :\n" << m << std::endl;

添加图片注释，不超过 140 字（可选）

inline Eigen::Block<Eigen::MatrixXf, -1, -1, false> Eigen::DenseBase<Eigen::MatrixXf>::block<int, int>(Eigen::Index startRow, Eigen::Index startCol, int blockRows, int blockCols)block<int, int>(Eigen::Index startRow, Eigen::Index startCol, int blockRows, int blockCols)

从矩阵中提取一个固定大小的子块。它接受起始行 (startRow)、起始列 (startCol)、子块的行数 (blockRows) 和列数 (blockCols) 作为参数

Eigen::Matrix2f b =
    m.block(1, 1, 2, 2);  // coping the middle part of matrix
std::cout << "Middle of 4x4 matrix :\n" << b << std::endl;

添加图片注释，不超过 140 字（可选）

将子块元素均变为0

m.block(1, 1, 2, 2) *= 0;  // change values in original matrix
std::cout << "Modified middle of 4x4 matrix :\n" << m << std::endl;

添加图片注释，不超过 140 字（可选）

第二行每个元素加3

 m.row(1).array() += 3;
std::cout << "Modified row of 4x4 matrix :\n" << m << std::endl;

添加图片注释，不超过 140 字（可选）

第三列每个元素除4

m.col(2).array() /= 4;
std::cout << "Modified col of 4x4 matrix :\n" << m << std::endl;

添加图片注释，不超过 140 字（可选）

Eigen的broadcasting操作

随机生成一个 2∗4 矩阵和一个 2∗1 的列向量

Eigen::MatrixXf mat = Eigen::MatrixXf::Random(2, 4);
std::cout << "Random 2x4 matrix :\n" << mat << std::endl;
Eigen::VectorXf v(2);  // column vector
v << 100, 100;

添加图片注释，不超过 140 字（可选）

把 v 的值加到矩阵的每一列对应的元素上

添加图片注释，不超过 140 字（可选）

线性回归基本理论

假设我们有一个数据集 {yi,xi1,…,xip}i=1n ，可以将 y 与 x 之间的线性关系表示为：

yi=β0⋅1+β1xi1+⋯+βpxip+ϵi=xiTβ+ϵi,i=1,…,n

其中， p 是自变量的维度， T 表示转置，因此 xiTβ 是向量 xi 与 β 的内积。

我们可以用矩阵的形式来重写该公式：

y=Xβ+ϵ

具体表示如下：

y=(y1y2⋮yn) X=(x1Tx2T⋮xnT)=(1x11…x1p1x21…x2p⋮⋮⋱⋮1xn1…xnp) β=(β0β1⋮βp) ϵ=(ϵ1ϵ2⋮ϵn)

当我们考虑简单线性回归时， p=1 ，此时公式如下：

yi=β0⋅1+β1xi+ϵi

线性回归任务的目标是找到满足前一个方程的参数向量。通常情况下，这样的线性方程组没有精确解，因此任务是在一些假设条件下估计满足这些方程的参数。最流行的估计方法之一是基于最小二乘原理的估计方法：最小化给定数据集中观测的因变量与线性函数预测的因变量之差的平方和。这被称为常规最小二乘( OLS )估计量。因此，任务可以用下面的公式来描述：

β^=arg⁡minβS(β)

目标函数S由下面的矩阵表示法给出：

S(β)=∑i=1n|yi−∑j=1pXijβj2|=||y−Xβ2||

在 X 矩阵的 p 列线性无关的情况下，这个极小化问题有唯一的解.我们可以通过求解法方程得到此解，具体如下：

β=(XTX)−1XTy

线性代数库可以用解析的方法直接求解此类方程，但它有一个显著的缺点-计算成本。在 y 和 x 维度较大的情况下，对计算机内存量和计算时间的要求太大，难以解决现实任务。

因此，通常用迭代方法来解决这个最小化任务。梯度下降(GD)就是这样一个算法的例子。即如果函数 S(β) 存在，并且在某点 β 的邻域内是可微的，那么在该点 β 沿着 S 的负梯度方向 S(β) 下降最快。

我们可以将我们的目标函数 S(β) 改为更适合迭代方法的形式。我们可以使用均方误差( MSE )函数，它衡量实际值和估测值之间的差异：

S(β)=1n∑i=1n(yi−Xiβ)2

在多元回归的情况下，我们对这个函数对每个 x 分量取偏导数 ∂S∂βj

所以，在线性回归的情况下，我们取如下的导数：

∂S∂β0=2n∑i=1n(Xiβ−yi)

∂S∂β1=2n∑i=1n(Xiβ−yi)Xi

整个算法的描述如下：

1. 用零初始化β。

2. 定义学习率参数的值，该参数控制我们在学习过程中调整参数的程度。

3. 计算以下β的值：

β0=β0−γ2n∑i=1n(Xiβ−yi)

β1=β1−γ2n∑i=1n(Xiβ−yi)Xi

4. 重复步骤1-3多次，直到均方误差（MSE）值达到合理的数值。

利用Eigen的正规方程和共轭梯度法实现

#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>

#include <algorithm>
#include <iostream>
#include <random>
#include <vector>

using Eigen::MatrixXf;

std::pair<Eigen::MatrixXf, Eigen::MatrixXf> GenerateData(size_t n) {
  std::vector<float> x_data(n);
  std::iota(x_data.begin(), x_data.end(), 0);
  std::vector<float> y_data(n);
  std::iota(y_data.begin(), y_data.end(), 0);

  // 添加随机噪声
  std::random_device rd;
  std::mt19937 re(rd());
  std::uniform_real_distribution<float> dist(-1.5f, 1.5f);

  for (auto& x : x_data) {
    x += dist(re);  // 添加噪声
  }

  for (auto& y : y_data) {
    y += dist(re);  // 添加噪声
  }

  // 生成数据集
  Eigen::Map<Eigen::MatrixXf> x(x_data.data(), static_cast<Eigen::Index>(n), 1);
  Eigen::Map<Eigen::MatrixXf> y(y_data.data(), static_cast<Eigen::Index>(n), 1);

  return {x, y};
}



int main() {
  size_t n = 1000;
  // 生成训练数据
  Eigen::MatrixXf x1, y;
  std::tie(x1, y) = GenerateData(n);
  Eigen::MatrixXf x0 = Eigen::MatrixXf::Ones(n, 1);
  // 对生成的y进行线性变换 y = b(4) + k(0.3)*x 使其具有线性关系
  y.array() *= 0.3f;
  y.array() += 4.f;
  Eigen::MatrixXf x(n, 2);
  x << x0, x1;

  // 使用共轭梯度法训练回归模型
  Eigen::LeastSquaresConjugateGradient<Eigen::MatrixXf> gd;
  gd.setMaxIterations(100);
  gd.setTolerance(0.001f);
  gdpute(x);
  Eigen::VectorXf b = gd.solve(y);
  std::cout << "Estimated parameters vector : \n" << b << std::endl;

  // 正规方程求解
  Eigen::VectorXf b_norm = (x.transpose() * x).ldlt().solve(x.transpose() * y);
  std::cout << "Estimated with normal equation parameters vector : " << b_norm
            << std::endl;

  // 预测
  Eigen::MatrixXf new_x(5, 2);
  new_x << 1, 1, 1, 2, 1, 3, 1, 4, 1, 5;
  auto new_y = new_x.array().rowwise() * b.transpose().array();
  std::cout << "Predicted values : \n" << new_y << std::endl;

  auto new_y_norm = new_x.array().rowwise() * b_norm.transpose().array();
  std::cout << "Predicted(norm) values : \n" << new_y_norm << std::endl;

  return 0;
};