Eigen:矩陣計算簡單用法(二)

查志強發表於2016-09-06
矩陣

【原文：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）矩陣的塊操作有兩種使用方法，其定義形式為：

  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


相關文章