Eigen C++開源矩陣計算工具——Eigen的簡單用法
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個模板引數,但是目前常用的只有前三個,如下所示:
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
struct traits<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
.......
其前三個引數分別表示矩陣元素的型別,行數和列數。typedef Matrix<double,Dynamic, Dynamic> MatrixXd; |
typedef Matrix< double , 3 , 1> Vector3d |
注意:
(1)Eigen中無論是矩陣還是陣列、向量,無論是靜態矩陣還是動態矩陣都提供預設建構函式,也就是你定義這些資料結構時都可以不用提供任何引數,其大小均由執行時來確定。
(2)矩陣的建構函式中只提供行列數、元素型別的構造引數,而不提供元素值的構造,對於比較小的、固定長度的向量提供初始化元素的定義,例如:
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:
#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:
#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開始)。可參看如下程式碼:
#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
Matrix3f m;
m << 1, 2, 3,
4, 5, 6,
7, 8, 9;
std::cout << m;
輸出結果為:
1 2 3 4 5 6 7 8 9程式碼段二(使用下標進行復制)
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) 使用“=”操作符操作動態矩陣時,如果左右邊的矩陣大小不等,則左邊的動態矩陣的大小會被修改為右邊的大小。例如下面的程式碼段:
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 3x36、如何選擇動態矩陣和靜態矩陣? Eigen對於這問題的答案是:對於小矩陣(一般大小小於16)的使用固定大小的靜態矩陣,它可以帶來比較高的效率,對於大矩陣(一般大小大於32)建議使用動態矩陣。
還需特別注意的是:如果特別大的矩陣使用了固定大小的靜態矩陣則可能造成棧溢位的問題
---------------------------------------------------------------------------------------------
本文主要是Eigen中矩陣和向量的算術運算,在Eigen中的這些算術運算過載了C++的+,-,*,所以使用起來非常方便。
1、矩陣的運算
Eigen提供+、-、一元操作符“-”、+=、-=,例如:
二元操作符+/-表示兩矩陣相加(矩陣中對應元素相加/減,返回一個臨時矩陣): B+C 或 B-C;
一元操作符-表示對矩陣取負(矩陣中對應元素取負,返回一個臨時矩陣): -C;
組合操作法+=或者-=表示(對應每隔元素都做相應操作):A += B 或者 A-=B
程式碼段1為矩陣的加減操作,程式碼如下:
#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也是一樣,例如下面的程式碼:
#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() 之類。
例如下面的程式碼所示:
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、矩陣相乘、矩陣向量相乘
矩陣的相乘,矩陣與向量的相乘也是使用操作符*,共有*和*=兩種操作符,其用法可以參考如下程式碼:
#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)矩陣的塊操作有兩種使用方法,其定義形式為:
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列的子矩陣,返回該子矩陣組成的臨時 矩陣物件,原矩陣的元素不變。
詳細使用情況,可參考下面的程式碼段:
#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開始。
下面的程式碼段用於演示獲取矩陣的指定行列:
#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 303)向量的塊操作,其實向量只是一個特殊的矩陣,但是Eigen也為它單獨提供了一些簡化的塊操作,如下三種形式:
獲取向量的前n個元素:vector.head(n);
獲取向量尾部的n個元素:vector.tail(n);
獲取從向量的第i個元素開始的n個元素:vector.segment(i,n);
其用法可參考如下程式碼段:
#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 | +++ | ++ |
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;
}
輸出結果為:
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
再分享一下我老師大神的人工智慧教程吧。零基礎!通俗易懂!風趣幽默!希望你也加入到我們人工智慧的隊伍中來!http://www.captainbed.net