C++中線性代數(shù)計算Eigen庫的使用教程詳解

更新時間：2023年12月30日 09:23:38 作者：yuanhao

Eigen是一個基于線性代數(shù)的C++模板庫,主要用于矩陣、向量、數(shù)值求解和相關(guān)算法,本文主要為大家簡單聊聊Eigen庫的使用,希望對大家有所幫助

前言

Eigen是一個基于線性代數(shù)的C++模板庫，主要用于矩陣、向量、數(shù)值求解和相關(guān)算法。Eigen庫有如下特點：

支持整數(shù)、浮點數(shù)、復(fù)數(shù)，使用模板編程，可以為特殊的數(shù)據(jù)結(jié)構(gòu)提供矩陣操作；
Open-CV自帶Eigen的接口；
支持逐元素、分塊和整體的矩陣操作。
支持使用Intel MKL加速部分功能。
支持多線程，對稀疏矩陣支持良好。
支持常用幾何運算，包括旋轉(zhuǎn)矩陣、矩陣變換等。

所以不論是做算法還是C++開發(fā)，使用好Eigen庫會讓你事半功倍。

由于本人是非專業(yè)的算法工程師，所以這里也只是做一些簡單入門介紹和使用。

類型定義

Eigen庫中的核心類就是Matrix，它表示矩陣，是一個模板類，定義如下：

template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
class Matrix {
    
}

其中前3個模板參數(shù)需要我們指定，后面3個參數(shù)使用默認(rèn)即可。

第一個參數(shù)typename _Scalar，表示該參數(shù)是一個數(shù)據(jù)類型，而不是像_Rows和_Cols這樣的變量。那么現(xiàn)在比如我想定義一個4x4的float矩陣，就可以如下定義：

Matrix<float,4,4> matrix44f;

這樣定義雖然簡單明了，但是有點長，所以在Eigen庫中提供了很多開源直接使用的模板類，比如：

Matrix4f，表示4x4的float矩陣；
Vector3f，表示列向量，長度3的float向量；
RowVector2i，表示行向量，長度2的int向量；

所以對于一些常見矩陣可以使用內(nèi)置的模板類。

除了定義已知大小的矩陣，Eigen庫還可以定義動態(tài)矩陣。那什么是動態(tài)矩陣和靜態(tài)矩陣呢？

靜態(tài)矩陣就是編譯時候就知道大小的矩陣，比如Matrix3d，在編譯時我們就知道這是一個3x3的double類型矩陣。

與之對應(yīng)的是動態(tài)矩陣，這種矩陣需要在運行后才知道大小，比如需要從vector<vector<double>>數(shù)據(jù)類型中構(gòu)建矩陣，只有當(dāng)代碼運行時才可以知道構(gòu)建的矩陣的行列個數(shù)。這時可以使用MatrixXd來表示任意大小的元素類型為double的矩陣變量，其中X就表示未知大小，即Dynamic動態(tài)的意思。

對于矩陣類型后面的后綴表示元素類型，Eigen庫定義了4種元素類型：

d：表示double類型；
f：表示float類型；
i：表示整數(shù)；
c：表示復(fù)數(shù)；

比如Matrix2c，表示一個2x2元素類型為復(fù)數(shù)的矩陣。

新建矩陣

對于默認(rèn)構(gòu)造函數(shù)的矩陣，分配了大小和內(nèi)存空間，但是沒有初始化矩陣元素。而對于基本數(shù)據(jù)類型，默認(rèn)初始化的話，其值是隨機(jī)的，不能使用。比如下面代碼：

    Eigen::Matrix2f m;
    std::cout << "m:" << std::endl << m << std::endl;
    Eigen::Vector4i n;
    std::cout << "n:" << std::endl << n << std::endl;
?
//運行結(jié)果：
m:
5.88604e-039            0
9.18341e-041  2.8026e-045
n:
-517829874
        -2
1952350234
   4199631

在Eigen中重載了std::cout的<<運算符，所以可以直接使用std::cout打印矩陣結(jié)果。所以定義完矩陣后，我們必須要對進(jìn)行初始化。對于矩陣的初始化，有如下幾種方式： 0. 直接賦值，可以使用<<進(jìn)行賦值。

    Eigen::Matrix2f m;  //2x2的float矩陣
    m << 1.9,2.884,     //進(jìn)行初始化，期望是4個值
         4.33,5.898;
    std::cout << "m:" << std::endl << m << std::endl;
?
    Eigen::Vector4i n;  //長度為4的列向量
    n << 3,4,2,6;
    std::cout << "n:" << std::endl << n << std::endl;
?
//運行結(jié)果：
m:
  1.9 2.884
 4.33 5.898
n:
3
4
2
6

這種方式只適合于矩陣的大小固定的情況，這種初始化方式很像使用列表初始化的數(shù)組，只有知道矩陣的大小，其構(gòu)造函數(shù)才能從左向右、從上向下來進(jìn)行初始化。比如可以定義動態(tài)大小的矩陣MatrixXd，這時就無法使用<<進(jìn)行直接賦值，比如：

    Eigen::MatrixXi xi;     //動態(tài)大小的int矩陣
    xi << 1,2,3,4,5,6;
    std::cout << "xi:" << std::endl << xi << std::endl;

因為這種定義的變量xi，并不知道其行列數(shù)目，可能是3x2、2x3等等，所以用6個int去賦值，構(gòu)造函數(shù)無法進(jìn)行構(gòu)造。

對于向量，可以在構(gòu)造的時候進(jìn)行初始化。

    Eigen::RowVector3d k(1,2,3);    //長度為3，類型為double的行向量
     std::cout << "k:" << std::endl << k << std::endl;

一些特殊函數(shù)，比如全部初始化為0，初始化為1，隨機(jī)數(shù)等。

Eigen::MatrixXi zero = MatrixXi::Zero(3,4);    //全部初始化為0的矩陣
std::cout << "------ zero ------" << std::endl << zero << std::endl;
Eigen::MatrixXi one = MatrixXi::Ones(4,4);      //全部初始化為1的矩陣
std::cout << "------ one ------" << std::endl << one << std::endl;
Eigen::MatrixXi i = MatrixXi::Identity(3,3);    //單位矩陣
std::cout << "------ i ------" << std::endl << i << std::endl;
Eigen::Matrix4f random = Matrix4f::Random();    //隨機(jī)矩陣
std::cout << "------ random ------" << std::endl << random << std::endl;

運行結(jié)果：

------ zero ------
0 0 0 0
0 0 0 0
0 0 0 0
------ one ------
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
------ i ------
1 0 0
0 1 0
0 0 1
------ random ------
-0.997497   0.170019    0.64568   0.421003
0.127171 -0.0402539    0.49321 0.0270699
-0.613392 -0.299417 -0.651784   -0.39201
0.617481   0.791925   0.717887 -0.970031

使用Eigen::Map將已有的內(nèi)存塊數(shù)據(jù)映射到矩陣，必須要指定需要映射的矩陣類型，是MatrixXd、VectorXd等等，在模板參數(shù)中指明。

看一下Eigen::Map函數(shù)的定義：

template< e> class Map
  : public MapBase<Map<PlainObjectType, MapOptions, StrideType> >
{
?
}

PlanObjectType：映射數(shù)據(jù)的等效矩陣類型。
MapOptions：指定指針對齊方式，以字節(jié)為單位，默認(rèn)為Unaligned，即未對齊。
StrideType：可選擇指定步幅。默認(rèn)情況下，Map采用內(nèi)存布局是一個普通的、連續(xù)的數(shù)組，可以通過指定步幅來定制化。

正常我們只需指明第一個參數(shù)即可，使用Map可以在無任何其他開銷的情況下，讓Eigen使用非Eigen數(shù)據(jù)結(jié)構(gòu)來對其進(jìn)行初始化，包括多種數(shù)據(jù)結(jié)構(gòu)，比如原生數(shù)組，比如double[]、float[]等，以及STL容器，比如std::vector、std::array等。

    std::vector<double> vec = {1,2,3,4,5,6};    //std::vector數(shù)據(jù)類型
    Eigen::Map<Eigen::VectorXd> vd(vec.data(),6);   //構(gòu)造成一個長度為6的列向量
    std::cout << "------ vd ------" << std::endl << vd << std::endl;
    Eigen::Map<Eigen::MatrixXd> xd(vec.data(),3,2); //構(gòu)造成一個3x2的矩陣
    std::cout << "------ xd ------" << std::endl << xd << std::endl;
    
    //運行結(jié)果：
------ vd ------
1
2
3
4
5
6
------ xd ------
1 4
2 5
3 6

在上面代碼中，我們使用vec.data()來獲取數(shù)vector中的數(shù)組數(shù)據(jù)，同時必須確定矩陣或者向量的大小。

對于矩陣來說，還有一種非常常見的初始化方式，就是一行一行或者一列一列進(jìn)行賦值。比如代碼：

Eigen::MatrixXd mat(3,2); //創(chuàng)建一個3x2的double類型矩陣
mat.col(0) << 1,2,3;    //對第一列進(jìn)行初始化
mat.col(1) << 4,5,6;    //對第二列進(jìn)行初始化

這種方式在構(gòu)建特殊矩陣時經(jīng)常用到。

矩陣索引

當(dāng)前矩陣的行數(shù)、列數(shù)和大小可以分別通過rows()、cols()和size()方法來獲取，在遍歷Eigen矩陣時最好獲取行數(shù)、列數(shù)來限制范圍。

同時索引下標(biāo)從0開始，矩陣元素的訪問可以通過matrix(i,j)來訪問第i行第j列的元素，對于向量還可以使用類似數(shù)組取值的vector[i]來訪問第i個元素。比如代碼：

    std::vector<double> vec = {1,2,3,4,5,6};
    Eigen::Map<Eigen::VectorXd> vd(vec.data(),6);   //創(chuàng)建一個長度為6的列向量
    std::cout << "------ vd ------" << std::endl << vd << std::endl;
    Eigen::Map<Eigen::MatrixXd> xd(vec.data(),3,2); //創(chuàng)建一個3x2的矩陣
    std::cout << "------ xd ------" << std::endl << xd << std::endl;
?
    for (int i = 0; i < xd.rows(); ++i) {
        for (int j = 0; j < xd.cols(); ++j) {
            std::cout << "xd.(" << i << "," << j << ")= " << xd(i,j) << std::endl;
        }
    }
?
    //運行結(jié)果：
------ xd ------
1 4
2 5
3 6
xd.(0,0)= 1
xd.(0,1)= 4
xd.(1,0)= 2
xd.(1,1)= 5
xd.(2,0)= 3
xd.(2,1)= 6

還可以利用block()函數(shù)來從Matrix中取出一個小矩陣來進(jìn)行處理，語法是matrix:block<p,q>(i,j)，表示從原矩陣的(i,j)位置開始，截取出pxq大小的子矩陣，比如測試代碼：

    Eigen::MatrixXi randomI = MatrixXi::Random(6,6);    //隨機(jī)值的6x6矩陣
    std::cout << "------ randomI ------" << std::endl << randomI << std::endl;
    Eigen::MatrixXi smallI = randomI.block<3,3>(1,1);   //從(1,1)處截取一個3x3的矩陣
    std::cout << "------ smallI ------" << std::endl << smallI << std::endl;
?
    //運行結(jié)果：
------ randomI ------
-13389 -12482   3334 -14515  12319  -1243
 -4442 -16231   3511   3528   7427  -8673
-11557 -16092 -10937   9283  14938  11869
-10948  -4002   5342   9915  13949  -9516
 16007   1037  -1613    651   1289   9163
 -1780   2332  -4846  -6490 -11720  11260
------ smallI ------
-16231   3511   3528
-16092 -10937   9283
 -4002   5342   9915

數(shù)學(xué)運算

Eigen庫重載了常用加減乘運算符，而對于矩陣的除法，我們一般是求逆然后再轉(zhuǎn)換為乘法，使用起來比較簡單，就不過多介紹了。

還有幾個常用的方法，必須要介紹一下，因為在項目中非常常見。

矩陣轉(zhuǎn)置。調(diào)用transpose()方法即可求出一個矩陣的轉(zhuǎn)置矩陣。
矩陣求逆。調(diào)用inverse()方法即可求出一個矩陣的逆矩陣。
共軛矩陣。調(diào)用conjugate()方法可以求出共軛矩陣。

求解線性最小二乘方程組

此外，Eigen庫的最主要一個作用可以直接求解最小二乘方程組，對于方程組Ax=b，如果沒有確定解，可以找到一個向量x，使得其誤差平方值最小。

我們可以回顧一下在線性回歸中，如何解出最佳擬合系數(shù)\theta，一般都是先構(gòu)建范德蒙矩陣，然后根據(jù)公式進(jìn)行計算。但是在Eigen中，可以直接使用bdcSvd類的solve()方法直接求解出最小二乘方程組的系數(shù)。

比如下面代碼中使用2種方法解最小二乘方程組：

    Eigen::MatrixXf X(4, 3);    //創(chuàng)建4x3的矩陣，表示有3個系數(shù)待求解，一共4個方程 
    Eigen::Vector4f Y;      //表示結(jié)果的長度為4的列向量
    X << 0, 0, 1,
    1, 1, 1,
    4, 2, 1,
    9, 3, 1;
    Y << 2, 0, 3, 7;
?
    std::cout << "------ X ------" << std::endl << X << std::endl;
    std::cout << "------ Y ------" << std::endl << Y << std::endl;
    //使用bdcSvd方法
    Eigen::MatrixXf abc = X.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
    std::cout << "------ abc ------" << std::endl << abc << std::endl;
    //使用公式法
    Eigen::MatrixXf abc1 = (X.transpose() * X).inverse() * (X.transpose()) * Y;
    std::cout << "------ abc1 ------" << std::endl << abc1 << std::endl;
?
    //運行結(jié)果：
------ X ------
0 0 1
1 1 1
4 2 1
9 3 1
------ Y ------
2
0
3
7
------ abc ------
 1.5
-2.7
 1.8
------ abc1 ------
     1.5
-2.70001
     1.8

到此這篇關(guān)于C++中線性代數(shù)計算Eigen庫的使用教程詳解的文章就介紹到這了,更多相關(guān)C++ Eigen庫內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: