1 矩陣基本運算簡介
Eigen重載了+,-,*運算符。同時提供了一些方法如dot(),cross()等。對於矩陣類的運算符重載只支持線性運算，比如matrix1*matrix2是矩陣相乘，當然必須要滿足矩陣乘法規則。對於向量和標量的加法(vector+scalar)這里並不支持，關於非線性運算這里暫不介紹。
2 加減運算
矩陣加減運算中必須要保證左右矩陣的行列對應相等。此外更重要的一點是，矩陣的類型也必須一致，這里的矩陣運算並不支持隱式的類型轉換。矩陣運算中重載的運算符有：
- 二元運算符+：a+b
- 二元運算符-：a-b
- 一元運算符-：-a
- 復合運算符+=：a+=b
- 復合運算符-=：a-=b
下面是使用示例：
#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
3 和標量的乘法和除法
對於矩陣和標量的乘法和除法也非常簡單，重載的操作符如下：
- 二元運算符：matrixscalar
- 二元運算符：scalarmatrix
- 二元運算符/：matrix/scalar
- 復合運算符=：matrix=scalar
- 復合運算符/=：matrix/=scalar
下面是使用示例：
#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
4 表達式計算優化原則
關於矩陣表達式的計算這里有一點需要說明。在Eigen中，算數運算操作(比如+)並不會立即對表達式兩端進行求值，而僅僅只是返回一個“表達式”，這個表達式對計算的結果的表現進行簡單的描述，而真正的計算會等到最后才會進行。一般來說會等到操作運算符=執行時進行計算。這個機制會極大的優化矩陣的計算性能。請看下面的這個例子：
VectorXf a(50),b(50),c(50),d(50);
...
a=3*b+4*c+5*d;
雖然上面這個表達式中有多個運算符，但並不會使用多個循環對每個運算符左右兩邊的矩陣進行求值。而是簡化為下面這一個循環。
for(int i = 0;i < 50;++i)
a[i] = 3*b[i] + 4*c[i] + 5*d[i];
所以我們並不用擔心較大或者復雜的運算表達式會降低我們的運算效率。它只會給Eigen提供更多的優化機會。
5 轉置(Transposition)和共軛(conjugation)
對於矩陣的轉置(transpose)
,共軛(conjugate)
以及伴隨矩陣(adjoint--conjugate transose)a*可以使用transpose(),conjugate(),'adjoint()'函數求得。
#include <iostream>
#include "Eigen\Dense"
using namespace Eigen;
int main()
{
MatrixXcf a = MatrixXcf::Random(2,2);
std::cout<<"Matrix a=\n"<<a<<std::endl;
std::cout<<"Here is the matrix a^T\n"<<a.transpose()<<std::endl;
std::cout<<"Here is the conjugate of a\n"<<a.conjugate()<<std::endl;
std::cout<<"Here is the matrix a^*\n"<<a.adjoint()<<std::endl;
}
執行結果如下：
Matrix a=
(0.127171,-0.997497) (-0.0402539,0.170019)
(0.617481,-0.613392) (0.791925,-0.299417)
Here is the matrix a^T
(0.127171,-0.997497) (0.617481,-0.613392)
(-0.0402539,0.170019) (0.791925,-0.299417)
Here is the conjugate of a
(0.127171,0.997497) (-0.0402539,-0.170019)
(0.617481,0.613392) (0.791925,0.299417)
Here is the matrix a^*
(0.127171,0.997497) (0.617481,0.613392)
(-0.0402539,-0.170019) (0.791925,0.299417)
對於實數矩陣，conjugate()不會有任何動作。所以adjoint()=transpose().
需要說明的是，作為基本運算符transpose()和adjoint()會簡單的返回一個沒有做任何轉換的代理對象(proxy object).如果使用b=a.transpose(),會使結果轉換和寫入b同時進行。然后這種轉換和寫入同時也會引起如下問題。如果我們執行a=a.transpose(),由於在轉換前就開始寫入了，所以並不會將轉換后的結果寫入a中。
#include <iostream>
#include "Eigen\Dense"
using namespace Eigen;
using namespace std;
int main()
{
Matrix2i a;
a<<1,2,3,4;
cout<<"Here is the matrix a:\n"<<endl;
a=a.transpose(); //!!! do NOT do this!!!
cout<<"and the result of the aliasing effect:\n"<<a<<endl;
}
官方給的例程說執行上面的程序會發現轉換后a的矩陣等於轉換前的，但是我測試的結果是程序直接出錯並停止運行。
如果我們想對a本身進行轉換可使用transposeInPlace()函數。同樣的如果求伴隨矩陣的話可使用adjoinInPlace()函數。
#include <iostream>
#include "Eigen\Dense"
using namespace Eigen;
using namespace std;
int main()
{
MatrixXf a(2,3);
a<<1,2,3,
4,5,6;
cout<<"Here is the initial matrix a:\n"<<a<<endl;
a.transposeInPlace();
cout<<"and after being transposed:\n"<<a<<endl;
a.adjointInPlace();
cout<<"and (a^T)^*=\n"<<a<<endl;
}
執行結果如下：
Here is the initial matrix a:
1 2 3
4 5 6
and after being transposed:
1 4
2 5
3 6
and (a^T)^*=
1 2 3
4 5 6
6 矩陣-矩陣以及矩陣-向量相乘
由於向量屬於特殊的矩陣，所以我們只需考慮矩陣的相乘即可。矩陣和矩陣相乘可以使用如下兩種運算符：
- 二元運算符：ab
- 復合運算符=：a=b
下面是使用示例：
#include <iostream>
#include "Eigen\Dense"
using namespace Eigen;
using namespace std;
int main()
{
Matrix2d mat;
mat<<1,2,
3,4;
Vector2d u(1,-1),v(2,0);
cout<<"Here is mat*mat:\n"<<mat*mat<<endl;
cout<<"Here is mat*u:\n"<<mat*u<<endl;
cout<<"Here is u^T*mat:\n"<<u.transpose()*mat<<endl;
cout<<"Here is u^T*v:\n"<<u.transpose()*v<<endl;
cout<<"Here is u*v^T:\n"<<u*v.transpose()<<endl;
cout<<"Let's multiply mat by itsef"<<endl;
mat = mat*mat;
cout<<"Now mat is mat:\n"<<mat<<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 itsef
Now mat is mat:
7 10
15 22
在矩陣的乘法中我們不用擔心a=a*a會引起上面的使用別名問題(aliasing issues)，因為這里會自動的引入一個中間變量因此a=a*a相當於temp=a*a;a=temp.如果你知道你的計算不會引起使用別名問題，那么你可以使用noalias()函數去避免使用這個中間變量以增加運算速度。如c.noalias()+=a*b;.關於使用別名的問題可以在官網aliasing中查看更多信息。
對於BLAS用戶可能擔心運算性能的問題，但正如我們前面所說的c.noalias-=2*a.adjoint()*b會完全的進行優化只觸發一個函數調用。
7 點乘和叉乘
對於內積和外積的計算可以使用dot()和cross()函數。當然點乘會得到一個1×1的矩陣(u.adjoint()*v)。
下面是一個使用示例：
#include <iostream>
#include "Eigen\Dense"
using namespace Eigen;
using namespace std;
int main()
{
Vector3d v(1,2,3);
Vector3d w(0,1,2);
cout<<"DOt product:"<<v.dot(w)<<endl;
double dp = v.adjoint()*w; //automatic conversion of the inner product to a scalar
cout<<"Dot product via a matrix product: "<<dp<<endl;
cout<<"Cross product:\n"<<v.cross(w)<<endl;
}
計算結果如下所示：
DOt product:8
Dot product via a matrix product: 8
Cross product:
1
-2
1
需要記住叉積只能對3維向量使用，而點積可以對任意維的向量使用。對於復數，點積會對第一個向量首先進行共軛然后在和第二個向量相乘。
8 其他一些基本運算
Eigen同樣提供了其他的函數對矩陣的所有元素進行操作，比如sum(對矩陣所有元素求和)，product（全部元素相乘），maximum（求最大值）和minimum（求最小值）。
下面是一個示例：
#include <iostream>
#include "Eigen\Dense"
using namespace Eigen;
using namespace std;
int main()
{
Matrix2d mat;
mat<<1,2,
3,4;
cout<<"Here is mat.sum():\t\t"<<mat.sum()<<endl;
cout<<"Here is mat.prd():\t\t"<<mat.prod()<<endl;
cout<<"Here is mat.mean():\t\t"<<mat.mean()<<endl;
cout<<"Here is mat.minCoeff():\t\t"<<mat.minCoeff()<<endl;
cout<<"Here is mat.maxCoeff():\t\t"<<mat.maxCoeff()<<endl;
cout<<"Here is mat.trace():\t\t"<<mat.trace()<<endl;
}
執行結果如下：
Here is mat.sum(): 10
Here is mat.prd(): 24
Here is mat.mean(): 2.5
Here is mat.minCoeff(): 1
Here is mat.maxCoeff(): 4
Here is mat.trace(): 5
對於求矩陣的跡除了使用trace()還可以使用高效的a.diagonal().sum()。
對於求minCoeff和maxCoeff除了能求出最大值以外還能獲取相關的索引下標。
#include <iostream>
#include <cstddef>
#include "Eigen\Dense"
using namespace Eigen;
using namespace std;
int main()
{
Matrix3f m = Matrix3f::Random();
ptrdiff_t i,j;
float minOfM = m.minCoeff(&i,&j);
cout<<"Here is the matrix m:\n"<<m<<endl;
cout<<"Its minimum coefficient ("<<minOfM
<<") is at position ("<<i<<","<<j<<")"<<endl<<endl;
RowVector4i v = RowVector4i::Random();
int maxOfv=v.maxCoeff(&i);
cout<<"Here is the vector v: "<<v<<endl;
cout<<"Its maximun coefficient ("<<maxOfv
<<") is at position "<<i<<endl;
}
執行結果如下：
Here is the matrix m:
-0.997497 0.617481 -0.299417
0.127171 0.170019 0.791925
-0.613392 -0.0402539 0.64568
Its minimum coefficient (-0.997497) is at position (0,0)
Here is the vector v: 8080 -10679 11761 6897
Its maximun coefficient (11761) is at position 2
9 有效性檢查
Eigen會進行操作的有效性檢測。如果可能的話在編譯階段就會給出相關的錯誤信息，雖然這些信息看起來又臭又長。但是Eigen會將重要的信息使用大寫加下划線的形式寫出。比如：
Matrix3f m;
Vector4f v;
v=m*v; //Compile-time error:YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
當然對於大多數情況，比如檢查動態數組大寫。編譯器並不能在編譯階段檢查出錯誤。Eigen在運行中會使用斷言(assertions)進行檢測。這意味着，如果我們使用了"debug mode",那么程序在檢測到一個非法操作后會被錯誤信息打斷。如果沒有使用斷言機制，那么程序可能會遇到災難性的錯誤。
MatrixXf m(3,3);
VectorXf v(4);
v = m*v; //Run-time assertion failure here:"invalid matrix product"