矩阵计算工具——Eigen的简单用法

最新推荐文章于 2022-08-15 20:02:01 发布

Jared_LAU

最新推荐文章于 2022-08-15 20:02:01 发布

阅读量3.1k

点赞数 1

分类专栏：编程语言与技巧文章标签： C++ 矩阵计算

编程语言与技巧专栏收录该内容

11 篇文章 0 订阅

订阅专栏

Eigen非常方便矩阵操作，当然它的功能不止如此，由于本人只用到了它的矩阵相关操作，所以这里只给出了它的一些矩阵相关的简单用法，以方便快速入门。矩阵操作在算法研究过程中，非常重要，例如在图像处理中二维高斯拟合求取光斑中心时使用Eigen提供的矩阵算法，差不多十来行代码即可实现，具体可见：http://blog.csdn.net/hjx_1000/article/details/8490653

Eigen的下载与安装，可参考下面两个博客：

http://blog.csdn.net/hjx_1000/article/details/8477522

或者：http://blog.csdn.net/abcjennifer/article/details/7781936；

Eigen帮助文档的地址：http://eigen.tuxfamily.org/dox/pages.html，本文中很多例子也是直接摘自这些帮助文档，

另外关于Eigen的论坛可以访问http://forum.kde.org/viewforum.php?f=74

Eigen用源码的方式提供给用户使用，在使用时只需要包含Eigen的头文件即可进行使用。

之所以采用这种方式，是因为Eigen采用模板方式实现，由于模板函数不支持分离编译，所以只能提供源码而不是动态库的方式供用户使用，不过这也也更方面用户使用和研究。关于模板的不支持分离编译的更多内容，请参考：http://blog.csdn.net/hjx_1000/article/details/8093701

1、矩阵的定义

Eigen中关于矩阵类的模板函数中，共有6个模板参数，但是目前常用的只有前三个，如下所示：

[cpp] view plain copy

template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct traits<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
.......

其前三个参数分别表示矩阵元素的类型，行数和列数。

矩阵定义时可以使用Dynamic来表示矩阵的行列数为未知，例如：

 
       typedef Matrix<double,Dynamic, Dynamic> MatrixXd; 
      

在Eigen中也提供了很多常见的简化定义形式，例如：

typedef Matrix< double , 3 , 1> Vector3d

注意：

（1）Eigen中无论是矩阵还是数组、向量，无论是静态矩阵还是动态矩阵都提供默认构造函数，也就是你定义这些数据结构时都可以不用提供任何参数，其大小均由运行时来确定。

（2）矩阵的构造函数中只提供行列数、元素类型的构造参数，而不提供元素值的构造，对于比较小的、固定长度的向量提供初始化元素的定义，例如：

[cpp] view plain copy

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

2、动态矩阵和静态矩阵

动态矩阵是指其大小在运行时确定，静态矩阵是指其大小在编译时确定，在Eigen中并未这样称呼矩阵。具体可见如下两段代码：

代码段1：

[cpp] view plain copy

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
MatrixXd m = MatrixXd::Random(3,3);
m = (m + MatrixXd::Constant(3,3,1.2)) * 50;
cout << "m =" << endl << m << endl;
VectorXd v(3);
v << 1, 2, 3;
cout << "m * v =" << endl << m * v << endl;
}

代码段2：

[cpp] view plain copy

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
Matrix3d m = Matrix3d::Random();
m = (m + Matrix3d::Constant(1.2)) * 50;
cout << "m =" << endl << m << endl;
Vector3d v(1,2,3);
cout << "m * v =" << endl << m * v << endl;
}

说明：

1）代码段1中MatrixXd表示任意大小的元素类型为double的矩阵变量，其大小只有在运行时被赋值之后才能知道； MatrixXd::Random(3,3)表示产生一个元素类型为double的3*3的临时矩阵对象。

2）代码段2中Matrix3d表示元素类型为double大小为3*3的矩阵变量，其大小在编译时就知道；

3）上例中向量的定义也是类似，不过这里的向量时列优先，在Eigen中行优先的矩阵会在其名字中包含有row，否则就是列优先。

4）向量只是一个特殊的矩阵，其一个维度为1而已，如：typedef Matrix< double , 3 , 1> Vector3d

3、矩阵元素的访问

在矩阵的访问中，行索引总是作为第一个参数，需注意Eigen中遵循大家的习惯让矩阵、数组、向量的下标都是从0开始。矩阵元素的访问可以通过（）操作符完成，例如m(2,3)即是获取矩阵m的第2行第3列元素（注意行列数从0开始）。可参看如下代码：

[cpp] view plain copy

#include <iostream>
#include <Eigen/Dense>
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 << "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;
}

其输出结果为：

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

针对向量还提供[]操作符，注意矩阵则不可如此使用，原因为：在C++中m[i, j]中逗号表达式 “i, j”的值始终都是“j”的值，即m[i, j]对于C++来讲就是m[j];

4、设置矩阵的元素

在Eigen中重载了"<<"操作符，通过该操作符即可以一个一个元素的进行赋值，也可以一块一块的赋值。另外也可以使用下标进行复制，例如下面两段代码：

代码段1

[cpp] view plain copy

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

输出结果为：

1 2 3
4 5 6
7 8 9

代码段二（使用下标进行复制）

[cpp] view plain copy

VectorXf m_Vector_A;
MatrixXf m_matrix_B;
int m_iN =-1;
bool InitData(int pSrc[100][100], int iWidth, int iHeight)
{
if (NULL == pSrc || iWidth <=0 || iHeight <= 0)
return false;
m_iN = iWidth*iHeight;
VectorXf tmp_A(m_iN);
MatrixXf tmp_B(m_iN, 5);
int i =0, j=0, iPos =0;
while(i<iWidth)
{
j=0;
while(j<iHeight)
{
tmp_A(iPos) = pSrc[i][j] * log((float)pSrc[i][j]);
tmp_B(iPos,0) = pSrc[i][j] ;
tmp_B(iPos,1) = pSrc[i][j] * i;
tmp_B(iPos,2) = pSrc[i][j] * j;
tmp_B(iPos,3) = pSrc[i][j] * i * i;
tmp_B(iPos,4) = pSrc[i][j] * j * j;
++iPos;
++j;
}
++i;
}
m_Vector_A = tmp_A;
m_matrix_B = tmp_B;
}

5、重置矩阵大小

当前矩阵的行数、列数、大小可以通过rows()，cols()和size()来获取，对于动态矩阵可以通过resize()函数来动态修改矩阵的大小.

需注意：

(1) 固定大小的矩阵是不能使用resize（）来修改矩阵的大小；

(2) resize（）函数会析构掉原来的数据，因此调用resize（）函数之后将不能保证元素的值不改变。
(3) 使用“=”操作符操作动态矩阵时，如果左右边的矩阵大小不等，则左边的动态矩阵的大小会被修改为右边的大小。例如下面的代码段：

[cpp] view plain copy

MatrixXf a(2,2);
std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;
MatrixXf b(3,3);
a = b;
std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;

输出结果为：

a is of size 2x2
a is now of size 3x3

6、如何选择动态矩阵和静态矩阵？

Eigen对于这问题的答案是：对于小矩阵（一般大小小于16）的使用固定大小的静态矩阵，它可以带来比较高的效率，对于大矩阵（一般大小大于32）建议使用动态矩阵。

还需特别注意的是：如果特别大的矩阵使用了固定大小的静态矩阵则可能造成栈溢出的问题

---------------------------------------------------------------------------------------------

本文主要是Eigen中矩阵和向量的算术运算，在Eigen中的这些算术运算重载了C++的+，-，*，所以使用起来非常方便。

1、矩阵的运算

Eigen提供+、-、一元操作符“-”、+=、-=，例如：

二元操作符+/-表示两矩阵相加（矩阵中对应元素相加/减，返回一个临时矩阵）： B+C 或 B-C；

一元操作符-表示对矩阵取负（矩阵中对应元素取负，返回一个临时矩阵）： -C;

组合操作法+=或者-=表示（对应每隔元素都做相应操作）：A += B 或者 A-=B

代码段1为矩阵的加减操作，代码如下：

[cpp] view plain copy

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
a << 1, 2,
3, 4;
MatrixXd b(2,2);
b << 2, 3,
1, 4;
std::cout << "a + b =\n" << a + b << std::endl;
std::cout << "a - b =\n" << a - b << std::endl;
std::cout << "Doing a += b;" << std::endl;
a += b;
std::cout << "Now a =\n" << a << std::endl;
Vector3d v(1,2,3);
Vector3d w(1,0,0);
std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}

输出结果为：

a + b =
3 5
4 8
a - b =
-1 -1
 2  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

另外，矩阵还提供与标量（单一个数字）的乘除操作，表示每个元素都与该标量进行乘除操作。例如：

二元操作符*在：A*a中表示矩阵A中的每隔元素都与数字a相乘，结果放在一个临时矩阵中，矩阵的值不会改变。

对于a*A、A/a、A*=a、A /=a也是一样，例如下面的代码：

[cpp] view plain copy

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
a << 1, 2,
3, 4;
Vector3d v(1,2,3);
std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;
std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;
std::cout << "Doing v *= 2;" << std::endl;
v *= 2;
std::cout << "Now v =\n" << v << std::endl;
}

输出结果为：

a * 2.5 =
2.5  5
7.5 10
0.1 * v =
0.1
0.2
0.3
Doing v *= 2;
Now v =
2
4
6

需要注意：

在Eigen中，算术操作例如 “操作符+”并不会自己执行计算操作，他们只是返回一个“算术表达式对象”，而实际的计算则会延迟到后面的赋值时才进行。这些不影响你的使用，它只是为了方便Eigen的优化。

2、求矩阵的转秩、共轭矩阵、伴随矩阵。

可以通过成员函数transpose(), conjugate(),和 adjoint()来完成，注意这些函数返回操作后的结果，而不会对原矩阵的元素进行直接操作，如果要让原矩阵的进行转换，则需要使用响应的InPlace函数，例如：transposeInPlace() 、 adjointInPlace() 之类。

例如下面的代码所示：

[cpp] view plain copy

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix a\n" << a << endl;
cout << "Here is the matrix a^T\n" << a.transpose() << endl;
cout << "Here is the conjugate of a\n" << a.conjugate() << endl;
cout << "Here is the matrix a^*\n" << a.adjoint() << endl;

输出结果为：

Here is the matrix a
 (-0.211,0.68) (-0.605,0.823)
 (0.597,0.566)  (0.536,-0.33)
Here is the matrix a^T
(-0.211,0.68) (0.597,0.566)
(-0.605,0.823) (0.536,-0.33)
Here is the conjugate of a
 (-0.211,-0.68) (-0.605,-0.823)
 (0.597,-0.566)    (0.536,0.33)
Here is the matrix a^*
(-0.211,-0.68) (0.597,-0.566)
(-0.605,-0.823)   (0.536,0.33)

3、矩阵相乘、矩阵向量相乘

矩阵的相乘，矩阵与向量的相乘也是使用操作符*，共有*和*=两种操作符，其用法可以参考如下代码：

[cpp] view plain copy

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d mat;
mat << 1, 2,
3, 4;
Vector2d u(-1,1), v(2,0);
std::cout << "Here is mat*mat:\n" << mat*mat << std::endl;
std::cout << "Here is mat*u:\n" << mat*u << std::endl;
std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl;
std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl;
std::cout << "Here is u*v^T:\n" << u*v.transpose() << std::endl;
std::cout << "Let's multiply mat by itself" << std::endl;
mat = mat*mat;
std::cout << "Now mat is mat:\n" << mat << std::endl;
}

输出结果为：

Here is mat*mat:
 7 10
15 22
Here is mat*u:
1
1
Here is u^T*mat:
2 2
Here is u^T*v:
-2
Here is u*v^T:
-2 -0
 2  0
Let's multiply mat by itself
Now mat is mat:
 7 10
15 22

--------------------------------------------------------------------------------------------

本节主要涉及Eigen的块操作以及QR分解，Eigen的QR分解非常绕人，搞了很久才搞明白是怎么回事，最后是一个使用Eigen的矩阵操作完成二维高斯拟合求取光点的代码例子，关于二维高斯拟合求取光点的详细内容可参考：http://blog.csdn.net/hjx_1000/article/details/8490653

1、矩阵的块操作

1）矩阵的块操作有两种使用方法，其定义形式为：

[cpp] view plain copy

matrix.block(i,j,p,q); （1）
matrix.block<p,q>(i,j); （2）

定义（1）表示返回从矩阵的(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象，原矩阵的元素不变。

定义（2）中block（p, q）可理解为一个p行q列的子矩阵，该定义表示从原矩阵中第(i, j)开始，获取一个p行q列的子矩阵，返回该子矩阵组成的临时矩阵对象，原矩阵的元素不变。

详细使用情况，可参考下面的代码段：

[cpp] view plain copy

#include <Eigen/Dense>
#include <iostream>
using namespace std;
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 << "Block in the middle" << endl;
cout << m.block<2,2>(1,1) << endl << endl;
for (int i = 1; i <= 3; ++i)
{
cout << "Block of size " << i << "x" << i << endl;
cout << m.block(0,0,i,i) << endl << endl;
}
}

输出的结果为：

Block in the middle
 6  7
10 11

Block of size 1x1
1

Block of size 2x2
1 2
5 6

Block of size 3x3
 1  2  3
 5  6  7
 9 10 11

通过上述方式获取的子矩阵即可以作为左值也可以作为右值，也就是即可以用这个子矩阵给其他矩阵赋值，也可以给这个子矩阵对象赋值。

2）矩阵也提供了获取其指定行/列的函数，其实获取某行/列也是一种特殊的获取子块。可以通过 .col()和 .row()来完成获取指定列/行的操作，参数为列/行的索引。
注意：
（1）需与获取矩阵的行数/列数的函数（ rows(), cols() ）的进行区别，不要弄混淆。
（2）函数参数为响应行/列的索引，需注意矩阵的行列均以0开始。
下面的代码段用于演示获取矩阵的指定行列：

[cpp] view plain copy

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::MatrixXf m(3,3);
m << 1,2,3,
4,5,6,
7,8,9;
cout << "Here is the matrix m:" << endl << m << endl;
cout << "2nd Row: " << m.row(1) << endl;
m.col(2) += 3 * m.col(0);
cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
cout << m << endl;
}

输出结果为：

Here is the matrix m:
1 2 3
4 5 6
7 8 9
2nd Row: 4 5 6
After adding 3 times the first column into the third column, the matrix m is:
 1  2  6
 4  5 18
 7  8 30

3）向量的块操作，其实向量只是一个特殊的矩阵，但是Eigen也为它单独提供了一些简化的块操作，如下三种形式：
获取向量的前n个元素：vector.head(n);
获取向量尾部的n个元素：vector.tail(n);
获取从向量的第i个元素开始的n个元素：vector.segment(i,n);
其用法可参考如下代码段：

[cpp] view plain copy

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::ArrayXf v(6);
v << 1, 2, 3, 4, 5, 6;
cout << "v.head(3) =" << endl << v.head(3) << endl << endl;
cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl;
v.segment(1,4) *= 2;
cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
}

输出结果为：

v.head(3) =
1
2
3

v.tail<3>() = 
4
5
6

after 'v.segment(1,4) *= 2', v =
1
4
6
8
10
6

2、QR分解
Eigen的QR分解非常绕人，它总共提供了下面这些矩阵的分解方式：

Decomposition	Method	Requirements on the matrix	Speed	Accuracy
PartialPivLU	partialPivLu()	Invertible	++	+
FullPivLU	fullPivLu()	None	-	+++
HouseholderQR	householderQr()	None	++	+
ColPivHouseholderQR	colPivHouseholderQr()	None	+	++
FullPivHouseholderQR	fullPivHouseholderQr()	None	-	+++
LLT	llt()	Positive definite	+++	+
LDLT	ldlt()	Positive or negative semidefinite	+++	++

由于我只用到了QR分解，而且Eigen的QR分解开始使用时确实不容易入手，因此这里只提供了householderQR的分解方式的演示代码：

[cpp] view plain copy

void QR2()
{
Matrix3d A;
A<<1,1,1,
2,-1,-1,
2,-4,5;
HouseholderQR<Matrix3d> qr;
qr.compute(A);
MatrixXd R = qr.matrixQR().triangularView<Upper>();
MatrixXd Q = qr.householderQ();
std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl;
std::cout << "A "<< std::endl <<A << std::endl << std::endl;
std::cout <<"qr.matrixQR()"<< std::endl << qr.matrixQR() << std::endl << std::endl;
std::cout << "R"<< std::endl <<R << std::endl << std::endl;
std::cout << "Q "<< std::endl <<Q << std::endl << std::endl;
std::cout <<"Q*R" << std::endl <<Q*R << std::endl << std::endl;
}