Eigen:矩阵计算简单用法(二)

【原文：http://blog.sina.com.cn/s/blog_691fc8920102v02q.html】

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

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

 

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

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

 

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

1、矩阵的运算

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

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

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

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

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

2. #include   
3. #include   
4. using namespace Eigen;  
5. int main()  
6. {  
7. Matrix2d a;  
8. a << 1, 2,  
9. 3, 4;  
10. MatrixXd b(2,2);  
11. b << 2, 3,  
12. 1, 4;  
13. std::cout << "a + b =\n" << a + b << std::endl;  
14. std::cout << "a - b =\n" << a - b << std::endl;  
15. std::cout << "Doing a += b;" << std::endl;  
16. a += b;  
17. std::cout << "Now a =\n" << a << std::endl;  
18. Vector3d v(1,2,3);  
19. Vector3d w(1,0,0);  
20. std::cout << "-v + w - v =\n" << -v + w - v << std::endl;  
21. } 

输出结果为：

 

 

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也是一样，例如下面的代码：

1. #include   
2. #include   
3. using namespace Eigen;  
4. int main()  
5. {  
6. Matrix2d a;  
7. a << 1, 2,  
8. 3, 4;  
9. Vector3d v(1,2,3);  
10. std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;  
11. std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;  
12. std::cout << "Doing v *= 2;" << std::endl;  
13. v *= 2;  
14. std::cout << "Now v =\n" << v << std::endl;  
15. }  

输出结果为：

 

 

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() 之类。

例如下面的代码所示：

1. MatrixXcf a = MatrixXcf::Random(2,2);  
2. cout << "Here is the matrix a\n" << a << endl;  
3. cout << "Here is the matrix a^T\n" << a.transpose() << endl;  
4. cout << "Here is the conjugate of a\n" << a.conjugate() << endl;  
5. 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、矩阵相乘、矩阵向量相乘

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

1. #include   
2. #include   
3. using namespace Eigen;  
4. int main()  
5. {  
6. Matrix2d mat;  
7. mat << 1, 2,  
8. 3, 4;  
9. Vector2d u(-1,1), v(2,0);  
10. std::cout << "Here is mat*mat:\n" << mat*mat << std::endl;  
11. std::cout << "Here is mat*u:\n" << mat*u << std::endl;  
12. std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl;  
13. std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl;  
14. std::cout << "Here is u*v^T:\n" << u*v.transpose() << std::endl;  
15. std::cout << "Let's multiply mat by itself" << std::endl;  
16. mat = mat*mat;  
17. std::cout << "Now mat is mat:\n" << mat << std::endl;  
18. }  

输出结果为：

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

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

   (i,j);    （2）  

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

 

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

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

1. #include   
2. #include   
3. using namespace std;  
4. int main()  
5. {  
6. Eigen::MatrixXf m(4,4);  
7. m << 1, 2, 3, 4,  
8. 5, 6, 7, 8,  
9. 9,10,11,12,  
10. 13,14,15,16;  
11. cout << "Block in the middle" << endl;  
12. cout << m.block<2,2>(1,1) << endl << endl;  
13. for (int i = 1; i <= 3; ++i)  
14. {  
15. cout << "Block of size " << i << "x" << i << endl;  
16. cout << m.block(0,0,i,i) << endl << endl;  
17. }  
18. }  

输出的结果为：

 

 

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开始。
下面的代码段用于演示获取矩阵的指定行列：

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

输出结果为：

 

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);
     其用法可参考如下代码段：

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

输出结果为：

 

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的分解方式的演示代码：

1. void QR2()  
2. {  
3.     Matrix3d A;  
4.     A<<1,1,1,  
5.         2,-1,-1,  
6.         2,-4,5;  
7.   
8.     HouseholderQR qr;  
9.     qr.compute(A);  
10.     MatrixXd R = qr.matrixQR().triangularView();  
11.     MatrixXd Q =  qr.householderQ();  
12.     std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl;  
13.     std::cout << "A "<< std::endl <<A << std::endl << std::endl;  
14.     std::cout <<"qr.matrixQR()"<< std::endl << qr.matrixQR() << std::endl << std::endl;  
15.     std::cout << "R"<< std::endl <<R << std::endl << std::endl;  
16.     std::cout << "Q "<< std::endl <<Q << std::endl << std::endl;  
17.     std::cout <<"Q*R" << std::endl <<Q*R << std::endl << std::endl;  
18. } 

输出结果为：

 

3、一个矩阵使用的例子：用矩阵操作完成二维高斯拟合，并求取光斑中心

下面的代码段是一个使用Eigen的矩阵操作完成二维高斯拟合求取光点的代码例子，关于二维高斯拟合求取光点的详细内容可参考：http://blog.csdn.net/hjx_1000/article/details/8490653

http://blog.csdn.net/houjixin/article/details/8490941

http://blog.csdn.net/houjixin/article/details/8492841

http://blog.csdn.net/houjixin/article/details/8494582