将从 http://eigen.tuxfamily.org/index.php?title=Main_Page#Download 下载下来的压缩包解压, 将其中的Eigen文件夹放入到开发的项目文件夹中, 在调用的时候只需要调用

#include "Eigen/Eigen” 或 #include "Eigen/Dense”, 其中不同的包叙述如下:

多种包

1. 矩阵的初始化

/*单个赋值法*/
int main()
{
  MatrixXd m(2,2);  //MatrixXd 代表 这个矩阵是double类型, X代表具体的行数都动态编译的
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << "Here is the matrix m:\n" << m << std::endl;
  VectorXd v(2);
  v(0) = 4;
  v(1) = v(0) - 1;
  std::cout << "Here is the vector v:\n" << v << std::endl;
  
  
  m.resize(4,1);      //MatrixXd还可以通过resize调整大小
}

/*统一赋值法*/
int main()
{
  VectorXd m(10); //定义一个向量
  m.setZero();
  m.setOnes();
  m.setConstant(value);+++
  m.setRandom();
  m.setLinSpaced(size, low, high);  //此处会将该向量大小改为size大小, 并令其在low到high范围内均匀取点, 即步长是(high-low)/size
  
  MatrixXd m(10,10);  //定义一个矩阵
  //下面的方法都会将m的大小改为rows, cols
  x.setZero(rows, cols);    
  x.setOnes(rows, cols);
  x.setConstant(rows, cols, value);
  x.setRandom(rows, cols);

  x.setIdentity();    //将m变为单位矩阵, 大小基于原矩阵大小
  x.setIdentity(rows, cols);       
}

/*基于Map方法, 将数组映射为向量, 需要注意的是其只支持一维数组或向量的初始化, 对于多维数组而言, 需要循环初始化*/

//连续映射
float data[] = {1,2,3,4};
Map<Vector3f> v1(data);       // uses v1 as a Vector3f object
Map<ArrayXf>  v2(data,3);     // uses v2 as a ArrayXf object
Map<Array22f> m1(data);       // uses m1 as a Array22f object
Map<MatrixXf> m2(data,2,2);   // uses m2 as a MatrixXf object

//间隔映射
float data[] = {1,2,3,4,5,6,7,8,9};
//其中InnerStride和OuterStride指的是步长对象, 而其泛型指的是步长的数值
Map<VectorXf,0,InnerStride<2> >  v1(data,3);                      // = [1,3,5], 此处的3代表v1是一个3*1维的行的向量
Map<VectorXf,0,InnerStride<> >   v2(data,3,InnerStride<>(3));     // = [1,4,7], 此处的第一个3指代的是v2的维度, 第二个3指代的是步长

//需要注意的是矩阵的维度不能超过按照步长最大可取的元素个数, 如果超过, 就会Undefined behavior
Map<VectorXf,0,InnerStride<2> >  v3(data,7);   
/*
结果为:
1
3
5
7
9
-1.97697e-22
-1.18826e+29
*/


/*
在用步长初始化矩阵时, 用的是OuterStride而不是向量的InnerStride

下面这两行矩阵初始化结果都为:
1 4 7
2 5 8
*/
Map<MatrixXf,0,OuterStride<3> >  m2(data,2,3);                    // 此处的2和3指代的m2是一个2*3维的矩阵
Map<MatrixXf,0,OuterStride<> >   m1(data,2,3,OuterStride<>(3));   // 也就是说步长即可在泛型中指定, 也可在初始化时传入一个步长对象


//对于多维矩阵, 我们可以通过
int row=10,col=5;
MatrixXf m(10,5);
float arr[]={1,2,3,4,5};
for(int i=0;i<row;i++ ){
  m.row(i)=Map<VectorXf>(arr,col);  
}
cout<<m<<endl;

//或者通过对多维vector的转换
vector<vector<float> > vec(5,vector<float>(5,1));
MatrixXf m(5,5)
for(int i=0;i<vec.size();i++){
  float *arr=vec[i].data();    //将一个float型数组的首地址指向vector中起始位置
  m.row(i)=Map<VectorXf>(arr,5); //不过此处需要注意的是不要对arr取sizeof, 因为此时sizeof失效了, 不再能得到真实的arr的大小了
}
cout<<m<<endl;


//最后是最常见的重载运算符赋值法
Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;
cout << m;

2. 矩阵的切片

/*
矩阵的切片在eigen中称作block, 可以通过如下方式
*/

int main()
{
  Eigen::MatrixXf m(4,4);
  m <<  1, 2, 3, 4,
        5, 6, 7, 8,
        9,10,11,12,
       13,14,15,16;
  
    cout<<m.block(1,0,2,1)<<endl; //从左到右有四个参数, 分别表示 (起点所在行, 起点所在列, 向下所取行数, 向右所取列数)
  /*
  结果是:
      5
      9
  即从第1行第0列开始, 向下取2行, 向右取1列, 即其本身那列, 换句话说就是起点为 (1,0) 的宽为2长为1的长方形
  */
  
  
  //除了上述block用法之外还可以用泛型block, 不过需要注意的是, 此处的泛型需要是显式定义的, 而不能传入一个参数, 即不能m.block<rows,cols>(1,0)这种形式, 这种形式的泛型会被忽略, 然后编译器会报错, 即找不到block这个方法
  cout<< m.block<2,3>(1,0) <<endl;    //即  block<向下所取行数, 向右所取列数>(起点所在行, 起点所在列)
  /*
  结果是:
      5  6  7
      9 10 11
  */
  
  
}

除却上面提到的外, 还有下面这些简易版的

快速矩阵初始化

topRows和bottomRows函数对于Vector也是可行的, 除此之外, 对于Vector, 有如下:

Vector版的快速矩阵初始化

3. 实现矩阵形式的最小二乘法

最小二乘法的矩阵形式为:
$Ax=b\rightarrow A^TAx=A^Tb\rightarrow x=(A^TA)^{-1}A^Tb$

$\color{red}{需要注意的是, 仅当 A 是满秩矩阵时, A^TA才可求逆, 否则其逆不存在。}$

那么第一步需要做的是对A进行QR分解, 即 $A=QR$ , 其中Q是正交矩阵, 即 $Q^{-1}=Q^T$ , 而R是右上三角矩阵, 即假如A是mn维, 则Q是 mm 维, R是 m*n 维, 只不过R只有右上角有值。即如下图所示。

QR分解

QR分解公式如下, 注意因为Q是正交矩阵所以 $Q^TQ=I$ :
$x = (A ^ { \top } A)^{-1}A ^ { \top } b \Rightarrow A ^ { \top } A x = A ^ { \top } b \Rightarrow R ^ { \top } R x = R ^ { \top } Q ^ { \top } b \Rightarrow R x = Q ^ { \top } b$

那么现在的问题就变成了 $x=R^{-1}Q^Tb$ , 那么我们可以更进一步, 将R进行LU分解, 也就是常见的高斯消去法, 在matlab中, 通过左除的形式, 即 R \ (Q'b) 这行代码。在python中, 通过 np.linalg.solve()函数。

LU分解如下图所示。 $A=LU$ 即分解为为L(下三角矩阵), 和U(上三角矩阵), 注意L中对角线上都是1。

LU分解

如果对于行数不等于列数的, 其LU分解为:

$\left[\begin{array} { r r r } { 0.68 } & { - 0.605 } & { - 0.0452 } \\ { - 0.211 } & { - 0.33 } & { 0.258 } \\ { 0.566 } & { 0.536 } & { - 0.27 } \\ { 0.597 } & { - 0.444 } & { 0.0268 } \\ { 0.823 } & { 0.108 } & { 0.904 } \end{array}\right]= \left[ \begin{array} { r r r r r } { 1 } & { 0 } & { 0 } & { 0 } & { 0 } \\ { - 0.299 } & { 1 } & { 0 } & { 0 } & { 0 } \\ { - 0.05 } & { 0.888 } & { 1 } & { 0 } & { 0 } \\ { 0.0296 } & { 0.705 } & { 0.768 } & { 1 } & { 0 } \\ { 0.285 } & { - 0.549 } & { 0.0436 } & { 0 } & { 1 } \end{array} \right]\left[\begin{array} { c c c } { 0.904 } & { 0.823 } & { 0.108 } \\ { 0 } & { 0.812 } & { 0.569 } \\ { 0 } & { 0 } & { - 1.1 } \\ { 0 } & { 0 } & { 0 } \\ { 0 } & { 0 } & { 0 } \end{array}\right]$

因此此时能得到L和U矩阵, 那么我们可以将R替换为 $R=L\cdot U$ , 那么对于 $Rx=Q^Tb$ , 我们就有:
$LUx=Q^Tb$
令 $Ux=y$ , 因此我们就可以先对 $Ly=Q^Tb$ 求解, 解出 $y$ 的值, 接着对 $Ux=y$ 求解, 即可得到 $x$ 的值。

不过需要注意的是, 一般而言, 对于非奇异矩阵(秩不是满秩), 譬如 $\left[ \begin{array} { l l } { 0 } & { 1 } \\ { 1 } & { 1 } \end{array} \right]$ , 是无法通过LU分解, 即高斯消去法来求得L和U矩阵的, 详细的说, 假设矩阵 $A$ 为(其中 $A_{11}$ 为k*k矩阵):
$A = \left[ \begin{array} { l l } { A _ { 11 } } & { A _ { 12 } } \\ { A _ { 21 } } & { A _ { 22 } } \end{array} \right]$
仅当矩阵 $A$ 满足如下定义时, 才可被LU分解:
$\operatorname { Rank } \left( A _ { 11 } \right) + k \geq \operatorname { Rank } \left( \left[ \begin{array} { l l } { A _ { 11 } } & { A _ { 12 } } \end{array} \right] \right) + \operatorname { Rank } \left( \left[ \begin{array} { c } { A _ { 11 } } \\ { A _ { 21 } } \end{array} \right] \right)$

所以说对于那些 $A^TA$ 不可求逆的最小二乘法, 即 $A$ 不是满秩的情况, 通过QR分解和LU分解的结合也是无法求出正确的解的。

对于这种情况, 我们只能通过梯度下降法来近似找到局部最优解。如果是矩阵形式的话, 可以通过 solve 函数, 无论是matlab中的solve函数, 还是python中的np.linalg.solve(), 还是c++ eigen库中的solve函数, 近似求解。

下面是eigen库通过先转换为QR分解再转换为LU分解的求解过程, 其中LU参考官方API

#include "Eigen/Eigen"
VectorXf main(const MatrixXf &A,const VectorXf &b){
    VectorXf x(b.rows());
    
    HouseholderQR<MatrixXf> QR=HouseholderQR<MatrixXf>(A);  //基于A矩阵创建QR对象
    MatrixXf Q=QR.householderQ(); //得到Q矩阵
    
    //因为A=QR, 而Q矩阵是正交矩阵, 所以R=Q'A, 然后用R矩阵构建LU对象
    FullPivLU<MatrixXf> LU(Q.transpose()*b);  //此时得到的LU矩阵对象相当于是L+U的混合, 需要将L和U从中分离开来
    

    //注意之前提到的对于非方阵矩阵的LU分解, 其中L矩阵的对角线上都是1, 所以可以认定为是一个单位矩阵和一个下三角矩阵的和
    MatrixXf L=MatrixXf::Identity(LU.matrixLU().rows(),LU.matrixLU().rows()); //设置一个单位矩阵
    L.block(0,0,LU.matrixLU().rows(),LU.matrixLU().cols()).triangularView<StrictlyLower>()=LU.matrixLU();  //加上下三角矩阵
    
    
    MatrixXf U=LU.matrixLU().triangularView<Upper>(); //而矩阵U对角线上不是1, 所以只需要取LU对象的上三角矩阵即可

    //此处可以输出看看L和U的维度
    cout<<L.rows()<<"-"<<L.cols()<<endl;  
    cout<<U.rows()<<"-"<<U.cols()<<endl;
    
    /*接下来我们首先需要解Ly=Q'b得到y, 然后再解 Ux=y, 得到x。 
    对于Ax=b这样的形式, 求解时的调用方式是, A.colPivHouseholderQr().solve(b)这样的形式, 返回的结果就是所要求的x。 因此我们先调用 y=L.colPivHouseholderQr().solve(Q'b), 再基于当前得到的结果调用 U.colPivHouseholderQr().solve(y), 两者合在一起即如下形式
    */
    x=U.colPivHouseholderQr().solve(L.colPivHouseholderQr().solve((Q.transpose()*b)));
    
    //此处可以输出看看x的维度
    cout<<x.rows()<<"-"<<x.cols()<<endl;
    
    //接下来我们来求一下loss
    float loss=(A*x-b).sum();
    
    cout<<loss<<endl;
    
    /*
    ###################################################################################
    那么如果发现矩阵A本身是非奇异矩阵, 即不是满秩矩阵该怎么办呢, 那就只能近似求解方程组了, 通过上面提到的.colPivHouseholderQr().solve()函数, 直接求解Ax=b中x的值, 即
    */
    x=A.colPivHouseholderQr().solve(b);
    
    return x;
}

上述解法只是其中一种解法, 无论是QR分解还是LU分解, eigen库中都有很多种不同的形式, 具体请参考eigen库的API, 归纳来说, 其总共提供了如下这几种矩阵分解方式。

矩阵分解方式

c++矩阵运算库eigen

c++矩阵运算库eigen

1. 矩阵的初始化

2. 矩阵的切片

3. 实现矩阵形式的最小二乘法