c++矩陣運算庫Eigen簡介

C++矩陣運算庫Eigen介紹

C++中的矩陣運算庫常用的有Armadillo,Eigen,OpenCV,ViennaCL,PETSc等。我自己在網上搜了一下不同運算庫的特點，最后選擇了Eigen。主要原因是Eigen體積較小，不用安裝也不用編譯，庫是以頭文件的形式給出，直接將它扔到我們自己的工程文件中即可，移植起來也無壓力。我們可以在Eigen官網下載源文件。

Eigen的HelloWorld

我這里使用的Eigen的版本為Eigen 3.3.3，源文件目錄如下：

可以直接用記事本打開INSTALL文件，里面有編譯和不編譯時分別怎么使用。我這里使用的方法是不對文件進行編譯。

1. 找個地方新建一個工程目錄，這里我在桌面上新建一個MatrixTest文件夾。
2. 將Eigen目錄中的Eigen文件夾拷貝到我們的MatrixTest目錄中。
3. 在MatrixTest中再建立一個main.cpp文件寫入如下代碼：

#include <iostream>
#include "Eigen\Core"

//import most common Eigen types
using namespace Eigen;

int main(int,char*[])
{
	Matrix3f m3;
	m3<<1,2,3,4,5,6,7,8,9;
	Matrix4f m4 = Matrix4f::Identity();
	Vector4i v4(1,2,3,4);
	std::cout<<"m3\n"<<m3<<"\nm4:\n"
	<<m4<<"\nv4:\n"<<v4<<std::endl;
}

4. 在MatrixTest目錄的地址欄中輸入cmd，然后用g++ main.cpp對文件進行編譯，再運行命令a。可以看到我們已經輸出了我們的矩陣了。

Eigen常規矩陣定義

Eigen的使用在官網上有詳細的介紹，這里對我學習過程中用到的基本操作進行介紹。首先是矩陣的定義。
在矩陣類的模板參數共有6個。一般情況下我們只需要關注前三個參數即可。前三個模板參數如下所示：

Matrix<typename Scalar,int RowsAtCompileTime,int ColsAtCompileTime>

1. Scalar參數為矩陣元素的類型，該參數可以是int,float,double等。
2. RowsAtCompileTime和ColsAtCompileTime是矩陣的行數和列數。

如Matrix<float,4,4> M44是定義一個4×4的矩陣，矩陣元素以float類型存儲。直接使用矩陣模板定義一個矩陣往往會覺得麻煩，Eigen提供了一些基本矩陣的別名定義，如typedef Matrix<float,4,4> Matrix4f.下面是一些內置的別名定義.來源於官方手冊：

typedef Matrix< std::complex<double> , 2 , 2 > Matrix2cd
typedef Matrix< std::complex<float> , 2 , 2 > Matrix2cf
typedef Matrix< double , 2 , 2 > Matrix2d
typedef Matrix< float , 2 , 2 > Matrix2f
typedef Matrix< int , 2 , 2 > Matrix2i
typedef Matrix< std::complex<double> , 3 , 3 > Matrix3cd
typedef Matrix< std::complex<float> , 3 , 3 > Matrix3cf
typedef Matrix< double , 3 , 3 > Matrix3d
typedef Matrix< float , 3 , 3 > Matrix3f
typedef Matrix< int , 3 , 3 > Matrix3i
typedef Matrix< std::complex<double> , 4 , 4 > Matrix4cd
typedef Matrix< std::complex<float> , 4 , 4 > Matrix4cf
typedef Matrix< double , 4 , 4 > Matrix4d
typedef Matrix< float , 4 , 4 > Matrix4f
typedef Matrix< int , 4 , 4 > Matrix4i
typedef Matrix< std::complex<double> , Dynamic , Dynamic > MatrixXcd
typedef Matrix< std::complex<float> , Dynamic , Dynamic > MatrixXcf
typedef Matrix< double , Dynamic , Dynamic > MatrixXd
typedef Matrix< float , Dynamic , Dynamic > MatrixXf
typedef Matrix< int , Dynamic , Dynamic > MatrixXi
typedef Matrix< std::complex<double> , 1, 2 > RowVector2cd
typedef Matrix< std::complex<float> , 1, 2 > RowVector2cf
typedef Matrix< double , 1, 2 > RowVector2d
typedef Matrix< float , 1, 2 > RowVector2f
typedef Matrix< int , 1, 2 > RowVector2i
typedef Matrix< std::complex<double> , 1, 3 > RowVector3cd
typedef Matrix< std::complex<float> , 1, 3 > RowVector3cf
typedef Matrix< double , 1, 3 > RowVector3d
typedef Matrix< float , 1, 3 > RowVector3f
typedef Matrix< int , 1, 3 > RowVector3i
typedef Matrix< std::complex<double> , 1, 4 > RowVector4cd
typedef Matrix< std::complex<float> , 1, 4 > RowVector4cf
typedef Matrix< double , 1, 4 > RowVector4d
typedef Matrix< float , 1, 4 > RowVector4f
typedef Matrix< int , 1, 4 > RowVector4i
typedef Matrix< std::complex<double> , 1, Dynamic > RowVectorXcd
typedef Matrix< std::complex<float> , 1, Dynamic > RowVectorXcf
typedef Matrix< double , 1, Dynamic > RowVectorXd
typedef Matrix< float , 1, Dynamic > RowVectorXf
typedef Matrix< int , 1, Dynamic > RowVectorXi
typedef Matrix< std::complex<double> , 2 , 1> Vector2cd
typedef Matrix< std::complex<float> , 2 , 1> Vector2cf
typedef Matrix< double , 2 , 1> Vector2d
typedef Matrix< float , 2 , 1> Vector2f
typedef Matrix< int , 2 , 1> Vector2i
typedef Matrix< std::complex<double> , 3 , 1> Vector3cd
typedef Matrix< std::complex<float> , 3 , 1> Vector3cf
typedef Matrix< double , 3 , 1> Vector3d
typedef Matrix< float , 3 , 1> Vector3f
typedef Matrix< int , 3 , 1> Vector3i
typedef Matrix< std::complex<double> , 4 , 1> Vector4cd
typedef Matrix< std::complex<float> , 4 , 1> Vector4cf
typedef Matrix< double , 4 , 1> Vector4d
typedef Matrix< float , 4 , 1> Vector4f
typedef Matrix< int , 4 , 1> Vector4i
typedef Matrix< std::complex<double> , Dynamic , 1> VectorXcd
typedef Matrix< std::complex<float> , Dynamic , 1> VectorXcf
typedef Matrix< double , Dynamic , 1> VectorXd
typedef Matrix< float , Dynamic , 1> VectorXf
typedef Matrix< int , Dynamic , 1> VectorXi

2 向量

向量被作為一種特殊的矩陣進行處理，即要么行為一要么列為一。除非顯式的說明為行向量，否則這里將向量默認為列向量。請看下面兩個別名定義：

typedef Matrix<float,3,1> Vector3f;
typedef Matrix<int,1,2> RowVector2i;

3 矩陣的動態空間分配

很多時候在程序的編譯階段也許我們並不知道矩陣具體的行列數，這時候使用動態控件分配就顯得很必要了。當我們給矩陣模板中參數RowsAtCompileTime或者ColsAtCompileTime參數指定為Dynamic時，表示該矩陣對應行或列為一個動態分配的值。下面是兩個動態矩陣的別名定義：

typedef Matrix<double,Dynamic,Dynamic> MatrixXd;
typedef Matrix<int,Dynamic,1> VectorXi;

4 矩陣的構建

經過上面的介紹以后，我們應該能定義一些基本的矩陣了。如：

Matrix3f a;   //定義一個float類型3×3固定矩陣a
MatrixXf b;   //定義一個float類型動態矩陣b(0×0)
Matrix<int,Dynamic,3> b;  //定義一個int類型動態矩陣(0×3)

對應動態矩陣，我們也可以在構造的時候給出矩陣所占用的空間，比如：

MatrixXf a(10,15);  //定義float類型10×15動態矩陣
VectorXf b(30); //定義float類型30×1動態矩陣(列向量)

為了保持一致性，我們也可以使用上面構造函數的形式定義一個固定的矩陣，即Matrix3f a(3,3)也是允許的。

上面矩陣在構造的過程中並沒有初始化，Eigen還為一些小的(列)向量提供了可以初始化的構造函數。如：

Vector2d a(5.0,6.0);
Vector3d b(5.0,6.0,7.0);
Vector4d c(5.0,6.0,7.0,8.0);

5 矩陣元素的訪問

Eigen提供了矩陣元素的訪問形式和matlab中矩陣的訪問形式非常相似，最大的不同是matlab中元素從1開始，而Eigen的矩陣中元素是從0開始訪問。對於矩陣，第一個參數為行索引，第二個參數為列索引。而對於向量只需要給出一個索引即可。

#include <iostream>
#include "Eigen\Core"

//import most common Eigen types
using namespace Eigen;

int main()
{
	MatrixXd m(2,2);
	m(0,0) = 3;
	m(1,0) = 2.5;
	m(0,1) = -1;
	m(1,1) = m(1,0) + m(0,1);
	
	std::cout<<"Hear 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;
}

輸出結果如下：

Hear is the matrix m:
  3  -1
2.5 1.5
Here is the vector v:
4
3

像m(index)這種訪問形式並不僅限於向量之中，對於矩陣也可以這樣訪問。這一點和matlab相同，我們知道在matlab中定義一個矩陣a(3,4)，如果我訪問a(5)相當於訪問a(2,2),這是因為在matlab中矩陣是按列存儲的。這里比較靈活，默認矩陣元素也是按列存儲的，但是我們也可以通過矩陣模板構造參數Options=RowMajor改變存儲方式(這個參數是我們還沒有提到的矩陣構造參數的第4個參數)。

6 一般初始化方法

對於矩陣的初始化，我們可以采用下面的方法方便且直觀的進行：

Matrix3f m;
m<<1,2,3,
   4,5,6,
   7,8,9;
std:cout<<m;

7 矩陣的大小

Eigen提供了rows(),cols(),size()方法來獲取矩陣的大小，同時也同了resize()方法從新改變動態數組的大小。

#include <iostream>
#include "Eigen\Core"

using namespace Eigen;

int main()
{
	MatrixXd m(2,5);
	m<<1,2,3,4,5,
	   6,7,8,9,10;
	m.resize(4,3);
	std::cout<<"The matrix m is:\n"<<m<<std::endl;
	std::cout<<"The matrix m is of size "
			 <<m.rows()<<"x"<<m.cols()<<std::endl;
	std::cout<<"It has "<<m.size()<<" coefficients"<<std::endl;
	VectorXd v(2);
	v<<1,2;
	v.resize(5);
	std::cout<<"The vector v is:\n"<<v<<std::endl;
	std::cout<<"The vector v is of size "<<v.size()<<std::endl;
	std::cout<<"As a matrix,v is of size "<<v.rows()
			 <<"x"<<v.cols()<<std::endl;
}

輸出結果如下：

The matrix m is:
           1            3            5
           6            8           10
           2            4 9.58787e-315
           7            9 1.17493e-309
The matrix m is of size 4x3
It has 12 coefficients
The vector v is:
           1
           2
1.17477e-309
 7.0868e-304
           0
The vector v is of size 5
As a matrix,v is of size 5x1

可以看到我們可以把矩陣任意的resize，但是resize后矩陣的元素會改變，如果resize后的矩陣比之前的大會出現一些未初始化的元素。如果被resize的矩陣按列存儲(默認)，那么resize命令和matlab中的reshape執行結果相同，只是matlab要求reshape的矩陣前后元素必須相同，也就是不允許resize后不能出現未初始化的元素。
對於固定大小的矩陣雖然也支持resize命令，但是resize后的大小只能是它本身的大小，否則就會報錯。因為resize前后如果矩陣大小一樣，就不會執行resize。如果我們不想在resize后改變矩陣的對應元素，那么可以使用conservativeResize()方法。對應上面程序中的m矩陣我們調用m.conservativeResize(4,3)后得到結果如下。其中因為行數增加了，增加的行會以未初始化的形式出現。

The matrix m is:
           1            2            3
           6            7            8
9.58787e-315   2.122e-314 1.52909e+249
           0            0 2.47039e+249

http://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html

8 賦值和大小變換

在Eigen中使用=可以直接將一個矩陣復制給另外一個矩陣，如果被復制的和賦值矩陣大小不一致，會自動對被復制矩陣執行resize函數。當然如果被復制的矩陣為固定矩陣當然就不會執行resize函數。當然也可以通過一些設置取消這個自動resize的過程。

using namespace Eigen;

int main()
{
	MatrixXf a(2,2);
	MatrixXf b(3,3);
	b<<1,2,3,
	   4,5,6,
	   7,8,9;
	a = b;
	std::cout<<a<<std::endl;
}

輸出結果：

1 2 3
4 5 6
7 8 9

9 固定矩陣和動態矩陣

什么時候使用固定矩陣什么時候使用動態矩陣呢？簡單的說：當矩陣尺寸比較小時使用固定矩陣（一般小於16），當矩陣尺寸較大時使用動態矩陣（一般大於32）。使用固定矩陣有更好的表現，它可以避免重復的動態內存分配，固定矩陣實際上是一個array。即Matrix4f mymatrix;事實上是float mymatrix[16];。所以這個是真的不花費任何運行時間。相反動態矩陣的建立需要在
heap中分配空間。即MatrixXf mymatrix(rows,colums);實際上是float *mymatrix = new float[rows*colums];.此外動態矩陣還要保存行和列的值。
當然固定矩陣也存在着顯而易見的弊端。當數組的大小過大時，固定數組的速度優勢就不那么明顯了，相反過大的固定數組有可能造成stack的溢出。這時候動態矩陣的靈活性就顯得十分重要了。

10 其他模板參數

最開始我們已經提到了建立一個矩陣一共有6個模板參數，其中有3個我們還沒有提到（其實第三個參數已經提到過了）。

Matrix<typename Scalar,
			 int RowsAtCompileTime,
			 int ColsAtCompileTime,
			 int Options=0,
			 int MaxRowsAtCompileTime = RowsAtCompileTime,
			 int MaxColsAtCompileTime = ColsAtCompileTime>

1. Options:這個參數決定了矩陣在存儲過程中實際是按行還是按列存儲。這個存儲方式在前面我們提到的矩陣變換時必須要注意。默認是按列存儲，我們可以顯示的使用Options=RowMajor讓矩陣實際按行存儲。如Matrix<float,2,3,RowMajor> a;.
2. MaxRowsAtCompileTime和MaxColsAtCompileTime:這兩個值是設定動態矩陣在分配空間時最大的行數和列數。如Matrix<float,Dynamic,Dynamic,0,3,4>;.

11 常規的矩陣typedef

我們前面給出了一些常用的矩陣typedef.其實可以總結如下：

1. MatrixNt對應的是Matrix<type,N,N>.比如MatrixXi對應的是Matrix<int,Dynamic,Dynamic>.
2. VectorNt對應的是Matrix<type,N,1>.比如Vector2f對應的是Matrix<float,2,1>.
3. RowVectorNt對應的是Matrix<type,1,N>.比如RowVector3d對應的是Matrix<double,1,3>.

其中:

1. N可以是2,3,4或者X(表示Dynamic).
2. t可以是i(int),f(float),d(double),cf(complex ),cd(complex ).只定義了這些類型的typedef並不表示只支持這些數據類型的運算。比如所有的整形類型的運算都支持（長的，短的，有符號的，無符號的）。