Eigen: C++開源矩陣計算庫

　　Eigen庫被分為一個Core模塊和幾個附加的模塊，每個模塊有一個相關的頭文件，使用該模塊時需要包含該頭文件，為了能便利的使用eigen的幾個模塊，Eigen提供了Dense和Eigen兩個頭文件，各個頭文件和模塊如下表

 

Module	Header file	Contents
Core	#include <Eigen/Core>	Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
Geometry	#include <Eigen/Geometry>	Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis)
LU	#include <Eigen/LU>	Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)
Cholesky	#include <Eigen/Cholesky>	LLT and LDLT Cholesky factorization with solver
Householder	#include <Eigen/Householder>	Householder transformations; this module is used by several linear algebra modules
SVD	#include <Eigen/SVD>	SVD decompositions with least-squares solver (JacobiSVD, BDCSVD)
QR	#include <Eigen/QR>	QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)
Eigenvalues	#include <Eigen/Eigenvalues>	Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver, ComplexEigenSolver)
Sparse	#include <Eigen/Sparse>	Sparse matrix storage and related basic linear algebra (SparseMatrix, SparseVector)  (see Quick reference guide for sparse matrices for details on sparse modules)
	#include <Eigen/Dense>	Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
	#include <Eigen/Eigen>	Includes Dense and Sparse header files (the whole Eigen library)

Eigen庫分為核心模塊和額外模塊兩個部分，

Eigen和Dense頭文件方便的同時包含了幾個頭文件可以供使用

--Core

有關矩陣和數組的類，由基本的線性代數（包含 三角形和自伴乘積相關），還有相應對數組的操作。

--Geometry

幾何學的類，有關轉換、平移、進位制、2d旋轉、3D旋轉（四元組和角軸相關）

--LU

邏輯單元的類，有關求逆，求行列式，LU分解解算器（FullPivLU,PartialPivLU）

--Cholesky

包含LLT和LDLT的Cholesky因式分解法（Cholesky分解是把一個對稱正定的矩陣表示成一個下三角矩陣L和其轉置的乘積的分解）

--Householder

Householder變換，這個模塊供幾何線性代數模塊使用

--SVD

奇異值分解，最小二乘解算器解決奇異值分解

--QR

QR分解求解，三種方法：HouseholderQR、ColPivHouseholderQR,FullPivhouseholderQR

--Eigenvalues

特征值和特征向量分解的方法：EigenSolver、SelfAdjointEigenSolver、ComplexEigenSolver

--Sparse

稀疏矩陣相關類，對於系數矩陣的存儲及相關線性代數

--Dense

包含：core、Geometry、LU、Cholesky、SVD、QR和Eigenvalues模塊（頭文件）

--Eigen

包含上述所有模塊（頭文件）

1、矩陣的定義

Eigen中關於矩陣類的模板函數中，共有6個參數，但是常用的只有前三個，如下所示：

template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
class Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >

......

前三個參數分別表示矩陣元素的類型，行數和列數。

矩陣定義時可以使用Dynamic 來表示矩陣的行列數是未知，例如：

typedef matrix<double, Dynamic, Dynamic> MarixXd;

在Eigen中也提供了很多常見的簡化定義形式，例如

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);

Vector2f v1(x, y);

Array3i   v2(x, y, z);
Vector4d  v3(x, y, z, w);
VectorXf  v5; // empty object ArrayXf v6(size);

//2D objects atrix4f m1; MatrixXf m5; // empty object MatrixXf m6(nb_rows, nb_columns);

   (3) Eigen矩陣的賦值方法

// 單點賦值
Eigen::Matrix<double, 3, 3> srcPoints;
    srcPoints(0, 0) = 4;
    srcPoints(1, 0) = 0;
    srcPoints(2, 0) = 0;
    srcPoints(0, 1) = 2;
//...

//<<運算符賦值
Vector3f  v1;     
v1 << x, y, z;

ArrayXf   v2(4);  
v2 << 1, 2, 3, 4;

//
Matrix3f  m1;   
m1 << 1, 2, 3,
      4, 5, 6,
      7, 8, 9;

int rows=5, cols=5;
MatrixXf m(rows,cols);
m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(),
     MatrixXf::Zero(3,cols-3),
     MatrixXf::Zero(rows-3,3),
     MatrixXf::Identity(rows-3,cols-3);
cout << m;

//對應的輸出： 1 2 3 0 0 4 5 6 0 0 7 8 9 0 0 0 0 0 1 0 0 0 0 0 1
（3）矩陣的預定義：

//vector x
x.setZero();
x.setOnes();
x.setConstant(value);
x.setRandom();
x.setLinSpaced(size, low, high);

VectorXf::Unit(4,1) == Vector4f(0,1,0,0)
                    == Vector4f::UnitY()

//matrix x
x.setZero(rows, cols);
x.setOnes(rows, cols);
x.setConstant(rows, cols, value);
x.setRandom(rows, cols);

x.setIdentity();//單位矩陣
x.setIdentity(rows, cols);     

（4）映射操作

　　可以將c語言類型的數組映射為矩陣或向量：

　　注意：

　　1.映射只適用於對一維數組進行操作，如果希望將二維數組映射為對應的矩陣，可以借助"mat.row(i)=Map<Vector> v(data[i],n)"進行循環來實現，其中data為二維數組名，n為數組第一維的維度。

　　2.應設置之后數組名和矩陣或向量名其實指向同一個地址，只是名字和用法不同，另外，對其中的任何一個進行修改都會修改另外一個

//連續映射
float data[] = {1,2,3,4};
Map<Vector3f> v1(data);       // uses v1 as a Vector3f object
Map<ArrayXf>  v2(data,3);     // uses v2 as a ArrayXf object
Map<Array22f> m1(data);       // uses m1 as a Array22f object
Map<MatrixXf> m2(data,2,2);   // uses m2 as a MatrixXf object

//間隔映射
float data[] = {1,2,3,4,5,6,7,8,9};
Map<VectorXf,0,InnerStride<2> >  v1(data,3);                      // = [1,3,5]
Map<VectorXf,0,InnerStride<> >   v2(data,3,InnerStride<>(3));     // = [1,4,7]
Map<MatrixXf,0,OuterStride<3> >  m2(data,2,3);                    // both lines     |1,4,7|
Map<MatrixXf,0,OuterStride<> >   m1(data,2,3,OuterStride<>(3));   // are equal to:  |2,5,8|

 

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;

//
vector.resize(size); vector.resizeLike(other_vector); vector.conservativeResize(size); // matrix.resize(nb_rows, nb_cols); matrix.resize(Eigen::NoChange, nb_cols); matrix.resize(nb_rows, Eigen::NoChange); matrix.resizeLike(other_matrix); matrix.conservativeResize(nb_rows, nb_cols);

輸出結果為：

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

 

6、如何選擇動態矩陣和靜態矩陣

Eigen對於這問題的答案是：對於小矩陣（一般大小小於16）的使用固定大小的靜態矩陣，它可以帶來比較高的效率，對於大矩陣（一般大小大於32）建議使用動態矩陣。

還需特別注意的是：如果特別大的矩陣使用了固定大小的靜態矩陣則可能造成棧溢出的問題.

 

7、矩陣的運算

Eigen重載了基礎的+ - * / += -= *= /= *可以表示標量和矩陣或者矩陣和矩陣

矩陣與矩陣的運算

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的優化

 

算術運算方法：

 

add  subtract	mat3 = mat1 + mat2;           mat3 += mat1; mat3 = mat1 - mat2;           mat3 -= mat1;
scalar product	mat3 = mat1 * s1;             mat3 *= s1;           mat3 = s1 * mat1; mat3 = mat1 / s1;             mat3 /= s1;
matrix/vector  products *	col2 = mat1 * col1; row2 = row1 * mat1;           row1 *= mat1; mat3 = mat1 * mat2;           mat3 *= mat1;
transposition  adjoint *	mat1 = mat2.transpose();      mat1.transposeInPlace(); mat1 = mat2.adjoint();        mat1.adjointInPlace();
dot product  inner product *	scalar = vec1.dot(vec2); scalar = col1.adjoint() * col2; scalar = (col1.adjoint() * col2).value();
outer product *	mat = col1 * col2.transpose();
norm  normalization *	scalar = vec1.norm();         scalar = vec1.squaredNorm() vec2 = vec1.normalized();     vec1.normalize(); // inplace
cross product *	#include <Eigen/Geometry> vec3 = vec1.cross(vec2);

Coefficient-wise & Array operators

Coefficient-wise operators for matrices and vectors:

Matrix API *	Via Array conversions
mat1.cwiseMin(mat2) mat1.cwiseMax(mat2) mat1.cwiseAbs2() mat1.cwiseAbs() mat1.cwiseSqrt() mat1.cwiseProduct(mat2) mat1.cwiseQuotient(mat2)	mat1.array().min(mat2.array()) mat1.array().max(mat2.array()) mat1.array().abs2() mat1.array().abs() mat1.array().sqrt() mat1.array() * mat2.array() mat1.array() / mat2.array()

Array operators:*

 

Arithmetic operators	array1 * array2     array1 / array2     array1 *= array2    array1 /= array2 array1 + scalar     array1 - scalar     array1 += scalar    array1 -= scalar
Comparisons	array1 < array2     array1 > array2     array1 < scalar     array1 > scalar array1 <= array2    array1 >= array2    array1 <= scalar    array1 >= scalar array1 == array2    array1 != array2    array1 == scalar    array1 != scalar
Trigo, power, and  misc functions  and the STL variants	array1.min(array2)            array1.max(array2)            array1.abs2() array1.abs()                  abs(array1) array1.sqrt()                 sqrt(array1) array1.log()                  log(array1) array1.exp()                  exp(array1) array1.pow(exponent)          pow(array1,exponent) array1.square() array1.cube() array1.inverse() array1.sin()                  sin(array1) array1.cos()                  cos(array1) array1.tan()                  tan(array1) array1.asin()                 asin(array1) array1.acos()                 acos(array1)

Reductions

Eigen provides several reduction methods such as: minCoeff() , maxCoeff() , sum() , prod() , trace() *, norm()*, squaredNorm() *, all() , and any() . All reduction operations can be done matrix-wise, column-wise or row-wise. Usage example:

5 3 1 mat = 2 7 8       9 4 6	mat.minCoeff();	1
	mat.colwise().minCoeff();	2 3 1
	mat.rowwise().minCoeff();	1 2 4

Typical use cases of all() and any():

//Typical use cases of all() and any(): if((array1 > 0).all()) ... // if all coefficients of array1 are greater than 0 ... if((array1 < array2).any()) ... // if there exist a pair i,j such that array1(i,j) < array2(i,j) ...

 

c++矩陣運算庫Eigen

　　Eigen 的簡介和下載安裝

最近因為要寫機械臂的運動學控制程序，因此難免就要在C++中進行矩陣運算。然而C++本身不支持矩陣運算，Qt庫中的矩陣運算能力也相當有限，因此考慮尋求矩陣運算庫Eigen的幫助。　

Eigen是一個C++線性運算的模板庫：他可以用來完成矩陣，向量，數值解等相關的運算。（Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.）

Eigen庫下載： http://eigen.tuxfamily.org/index.php?title=Main_Page

Eigen庫的使用相當方便，將壓縮包中的Eigen文件夾拷貝到項目目錄下，直接包含其中的頭文件即可使用，省去了使用Cmake進行編譯的煩惱。

Eigen官網有快速入門的參考文檔：http://eigen.tuxfamily.org/dox/group__QuickRefPage.html

 

　　Eigen簡單上手使用

要實現相應的功能只需要包含頭相應的頭文件即可：

Core	#include <Eigen/Core>	Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
Geometry	#include <Eigen/Geometry>	Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis)
LU	#include <Eigen/LU>	Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)
Cholesky	#include <Eigen/Cholesky>	LLT and LDLT Cholesky factorization with solver
Householder	#include <Eigen/Householder>	Householder transformations; this module is used by several linear algebra modules
SVD	#include <Eigen/SVD>	SVD decomposition with least-squares solver (JacobiSVD)
QR	#include <Eigen/QR>	QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)
Eigenvalues	#include <Eigen/Eigenvalues>	Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver,ComplexEigenSolver)
Sparse	#include <Eigen/Sparse>	Sparse matrix storage and related basic linear algebra (SparseMatrix, DynamicSparseMatrix,SparseVector)
	#include <Eigen/Dense>	Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
	#include <Eigen/Eigen>	Includes Dense and Sparse header files (the whole Eigen library)

 基本的矩陣運算只需要包含Dense即可

基本的數據類型：Array, matrix and vector

聲明：

#include<Eigen/Dense>

...
//1D objects Vector4d v4; Vector2f v1(x, y); Array3i v2(x, y, z); Vector4d v3(x, y, z, w); VectorXf v5; // empty object ArrayXf v6(size); //2D objects atrix4f m1; MatrixXf m5; // empty object MatrixXf m6(nb_rows, nb_columns);

賦值：

//    
Vector3f  v1;     v1 << x, y, z; ArrayXf v2(4); v2 << 1, 2, 3, 4; // Matrix3f m1; m1 << 1, 2, 3, 4, 5, 6, 7, 8, 9;

int rows=5, cols=5; MatrixXf m(rows,cols); m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(), MatrixXf::Zero(3,cols-3), MatrixXf::Zero(rows-3,3), MatrixXf::Identity(rows-3,cols-3); cout << m; //對應的輸出： 1 2 3 0 0 4 5 6 0 0 7 8 9 0 0 0 0 0 1 0 0 0 0 0 1

Resize矩陣：

//
vector.resize(size); vector.resizeLike(other_vector); vector.conservativeResize(size); // matrix.resize(nb_rows, nb_cols); matrix.resize(Eigen::NoChange, nb_cols); matrix.resize(nb_rows, Eigen::NoChange); matrix.resizeLike(other_matrix); matrix.conservativeResize(nb_rows, nb_cols);

元素讀取：

vector(i)；    
vector[i]；

matrix(i,j)；

矩陣的預定義：

//vector x x.setZero(); x.setOnes(); x.setConstant(value); x.setRandom(); x.setLinSpaced(size, low, high); VectorXf::Unit(4,1) == Vector4f(0,1,0,0) == Vector4f::UnitY() //matrix x x.setZero(rows, cols); x.setOnes(rows, cols); x.setConstant(rows, cols, value); x.setRandom(rows, cols); x.setIdentity();//單位矩陣 x.setIdentity(rows, cols); 

映射操作，可以將c語言類型的數組映射為矩陣或向量：

（注意：

1.映射只適用於對一維數組進行操作，如果希望將二維數組映射為對應的矩陣，可以借助"mat.row(i)=Map<Vector> v(data[i],n)"進行循環來實現，其中data為二維數組名，n為數組第一維的維度。

2.應設置之后數組名和矩陣或向量名其實指向同一個地址，只是名字和用法不同，另外，對其中的任何一個進行修改都會修改另外一個）

//連續映射 float data[] = {1,2,3,4}; Map<Vector3f> v1(data); // uses v1 as a Vector3f object Map<ArrayXf> v2(data,3); // uses v2 as a ArrayXf object Map<Array22f> m1(data); // uses m1 as a Array22f object Map<MatrixXf> m2(data,2,2); // uses m2 as a MatrixXf object //間隔映射 float data[] = {1,2,3,4,5,6,7,8,9}; Map<VectorXf,0,InnerStride<2> > v1(data,3); // = [1,3,5] Map<VectorXf,0,InnerStride<> > v2(data,3,InnerStride<>(3)); // = [1,4,7] Map<MatrixXf,0,OuterStride<3> > m2(data,2,3); // both lines |1,4,7| Map<MatrixXf,0,OuterStride<> > m1(data,2,3,OuterStride<>(3)); // are equal to: |2,5,8|

算術運算：

 

add  subtract	mat3 = mat1 + mat2;           mat3 += mat1; mat3 = mat1 - mat2;           mat3 -= mat1;
scalar product	mat3 = mat1 * s1;             mat3 *= s1;           mat3 = s1 * mat1; mat3 = mat1 / s1;             mat3 /= s1;
matrix/vector  products *	col2 = mat1 * col1; row2 = row1 * mat1;           row1 *= mat1; mat3 = mat1 * mat2;           mat3 *= mat1;
transposition  adjoint *	mat1 = mat2.transpose();      mat1.transposeInPlace(); mat1 = mat2.adjoint();        mat1.adjointInPlace();
dot product  inner product *	scalar = vec1.dot(vec2); scalar = col1.adjoint() * col2; scalar = (col1.adjoint() * col2).value();
outer product *	mat = col1 * col2.transpose();
norm  normalization *	scalar = vec1.norm();         scalar = vec1.squaredNorm() vec2 = vec1.normalized();     vec1.normalize(); // inplace
cross product *	#include <Eigen/Geometry> vec3 = vec1.cross(vec2);

Coefficient-wise & Array operators

Coefficient-wise operators for matrices and vectors:

Matrix API *	Via Array conversions
mat1.cwiseMin(mat2) mat1.cwiseMax(mat2) mat1.cwiseAbs2() mat1.cwiseAbs() mat1.cwiseSqrt() mat1.cwiseProduct(mat2) mat1.cwiseQuotient(mat2)	mat1.array().min(mat2.array()) mat1.array().max(mat2.array()) mat1.array().abs2() mat1.array().abs() mat1.array().sqrt() mat1.array() * mat2.array() mat1.array() / mat2.array()

Array operators:*

 

Arithmetic operators	array1 * array2     array1 / array2     array1 *= array2    array1 /= array2 array1 + scalar     array1 - scalar     array1 += scalar    array1 -= scalar
Comparisons	array1 < array2     array1 > array2     array1 < scalar     array1 > scalar array1 <= array2    array1 >= array2    array1 <= scalar    array1 >= scalar array1 == array2    array1 != array2    array1 == scalar    array1 != scalar
Trigo, power, and  misc functions  and the STL variants	array1.min(array2)            array1.max(array2)            array1.abs2() array1.abs()                  abs(array1) array1.sqrt()                 sqrt(array1) array1.log()                  log(array1) array1.exp()                  exp(array1) array1.pow(exponent)          pow(array1,exponent) array1.square() array1.cube() array1.inverse() array1.sin()                  sin(array1) array1.cos()                  cos(array1) array1.tan()                  tan(array1) array1.asin()                 asin(array1) array1.acos()                 acos(array1)

Reductions

Eigen provides several reduction methods such as: minCoeff() , maxCoeff() , sum() , prod() , trace() *, norm()*, squaredNorm() *, all() , and any() . All reduction operations can be done matrix-wise, column-wise or row-wise. Usage example:

5 3 1 mat = 2 7 8       9 4 6	mat.minCoeff();	1
	mat.colwise().minCoeff();	2 3 1
	mat.rowwise().minCoeff();	1 2 4

Typical use cases of all() and any():

//Typical use cases of all() and any(): if((array1 > 0).all()) ... // if all coefficients of array1 are greater than 0 ... if((array1 < array2).any()) ... // if there exist a pair i,j such that array1(i,j) < array2(i,j) ...

 

Sub-matrices

Read-write access to a column or a row of a matrix (or array):

mat1.row(i) = mat2.col(j);

mat1.col(j1).swap(mat1.col(j2));

Read-write access to sub-vectors:

Default versions	Optimized versions when the size  is known at compile time
vec1.head(n)	vec1.head<n>()	the first n coeffs
vec1.tail(n)	vec1.tail<n>()	the last n coeffs
vec1.segment(pos,n)	vec1.segment<n>(pos)	the n coeffs in the  range [pos : pos + n - 1]
Read-write access to sub-matrices:
mat1.block(i,j,rows,cols) (more)	mat1.block<rows,cols>(i,j) (more)	the rows x cols sub-matrix  starting from position (i,j)
mat1.topLeftCorner(rows,cols) mat1.topRightCorner(rows,cols) mat1.bottomLeftCorner(rows,cols) mat1.bottomRightCorner(rows,cols)	mat1.topLeftCorner<rows,cols>() mat1.topRightCorner<rows,cols>() mat1.bottomLeftCorner<rows,cols>() mat1.bottomRightCorner<rows,cols>()	the rows x cols sub-matrix  taken in one of the four corners
mat1.topRows(rows) mat1.bottomRows(rows) mat1.leftCols(cols) mat1.rightCols(cols)	mat1.topRows<rows>() mat1.bottomRows<rows>() mat1.leftCols<cols>() mat1.rightCols<cols>()	specialized versions of block()  when the block fit two corners

top

Miscellaneous operations

Reverse

Vectors, rows, and/or columns of a matrix can be reversed (see DenseBase::reverse(), DenseBase::reverseInPlace(),VectorwiseOp::reverse()).

vec.reverse()           mat.colwise().reverse()   mat.rowwise().reverse()

vec.reverseInPlace()

Replicate

Vectors, matrices, rows, and/or columns can be replicated in any direction (see DenseBase::replicate(),VectorwiseOp::replicate())

vec.replicate(times)                                          vec.replicate<Times>

mat.replicate(vertical_times, horizontal_times)               mat.replicate<VerticalTimes, HorizontalTimes>()

mat.colwise().replicate(vertical_times, horizontal_times)     mat.colwise().replicate<VerticalTimes, HorizontalTimes>()

mat.rowwise().replicate(vertical_times, horizontal_times)     mat.rowwise().replicate<VerticalTimes, HorizontalTimes>()

top

Diagonal, Triangular, and Self-adjoint matrices

(matrix world *)

Diagonal matrices

Operation	Code
view a vector as a diagonal matrix	mat1 = vec1.asDiagonal();
Declare a diagonal matrix	DiagonalMatrix<Scalar,SizeAtCompileTime> diag1(size); diag1.diagonal() = vector;
Access the diagonal and super/sub diagonals of a matrix as a vector (read/write)	vec1 = mat1.diagonal();        mat1.diagonal() = vec1;      // main diagonal vec1 = mat1.diagonal(+n);      mat1.diagonal(+n) = vec1;    // n-th super diagonal vec1 = mat1.diagonal(-n);      mat1.diagonal(-n) = vec1;    // n-th sub diagonal vec1 = mat1.diagonal<1>();     mat1.diagonal<1>() = vec1;   // first super diagonal vec1 = mat1.diagonal<-2>();    mat1.diagonal<-2>() = vec1;  // second sub diagonal
Optimized products and inverse	mat3  = scalar * diag1 * mat1; mat3 += scalar * mat1 * vec1.asDiagonal(); mat3 = vec1.asDiagonal().inverse() * mat1 mat3 = mat1 * diag1.inverse()

Triangular views

TriangularView gives a view on a triangular part of a dense matrix and allows to perform optimized operations on it. The opposite triangular part is never referenced and can be used to store other information.

Note

The .triangularView() template member function requires the  template keyword if it is used on an object of a type that depends on a template parameter; see The template and typename keywords in C++ for details.

Operation	Code
Reference to a triangular with optional  unit or null diagonal (read/write):	m.triangularView<Xxx>() Xxx = Upper, Lower, StrictlyUpper, StrictlyLower, UnitUpper,UnitLower
Writing to a specific triangular part: (only the referenced triangular part is evaluated)	m1.triangularView< Eigen::Lower>() = m2 + m3
Conversion to a dense matrix setting the opposite triangular part to zero:	m2 = m1.triangularView< Eigen::UnitUpper>()
Products:	m3 += s1 * m1.adjoint().triangularView< Eigen::UnitUpper>() * m2 m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView< Eigen::Lower>()
Solving linear equations:	L1.triangularView< Eigen::UnitLower>().solveInPlace(M2) L1.triangularView<Eigen:: Lower>().adjoint().solveInPlace(M3) U1.triangularView<Eigen:: Upper>().solveInPlace< OnTheRight>(M4)

 

8、求矩陣的轉置、共軛矩陣、伴隨矩陣、行列式、逆矩陣。

可以通過 成員函數transpose(), conjugate(),determinant(), adjoint()和inverse()來完成，注意這些函數返回操作后的結果，而不會對原矩陣的元素進行直接操作，如果要讓原矩陣的進行轉換，則需要使用響應的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;
cout << "determinant： " << c.determinant() << endl;
cout << "inverse:" << c.inverse() << 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)
......

矩陣的特征值: m.eigenvalues();

特征向量: m.eigenvector();

 

9、矩陣相乘、矩陣向量相乘

矩陣的相乘，矩陣與向量的相乘也是使用操作符*，共有*和*=兩種操作符，其用法可以參考如下代碼：

#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

10、點積和叉積

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
    //點積、叉積(針對向量的)
    Vector3d v(1,2,3);
    Vector3d w(0,1,2);
    std::cout<<v.dot(w)<<std::endl<<std::endl;
    std::cout<<w.cross(v)<<std::endl<<std::endl;
}
*/

 

11、矩陣的塊操作

   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 30

3）向量的塊操作，其實向量只是一個特殊的矩陣，但是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

12、計算特征值和特征向量

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    //特征向量、特征值
    std::cout << "Here is the matrix A:\n" << a << std::endl;
    SelfAdjointEigenSolver<Matrix2d> eigensolver(a);
    if (eigensolver.info() != Success) abort();
     std::cout << "特征值:\n" << eigensolver.eigenvalues() << std::endl;
     std::cout << "Here's a matrix whose columns are eigenvectors of A \n"
     << "corresponding to these eigenvalues:\n"
     << eigensolver.eigenvectors() << std::endl;
}

13、矩陣分解

矩陣分解 (decomposition, factorization)是將矩陣拆解為數個矩陣的乘積，可分為三角分解、滿秩分解、QR分解、Jordan分解和SVD（奇異值）分解等，常見的有三種：1)三角分解法 (Triangular Factorization)，2)QR 分解法 (QR Factorization)，3)奇異值分解法 (Singular Value Decompostion)。

Eigen總共提供了下面這些矩陣的分解方式：

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分解

   QR分解法是將矩陣分解成一個正規正交矩陣與上三角形矩陣,所以稱為QR分解法,與此正規正交矩陣的通用符號Q有關。

householderQR的分解方式的演示代碼：

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;
}

輸出為：

QR2(): HouseholderQR---------------------------------------------
A 
 1  1  1
 2 -1 -1
 2 -4  5

qr.matrixQR()
 -3   3  -3
0.5  -3   3
0.5  -1  -3

R
-3  3 -3
 0 -3  3
 0  0 -3

Q 
-0.333333 -0.666667 -0.666667
-0.666667 -0.333333  0.666667
-0.666667  0.666667 -0.333333

Q*R
 1  1  1
 2 -1 -1
 2 -4  5

矩陣的LU三角分解 

 

矩陣的SVD（奇異值分解）分解

 

奇異值分解 (singular value decomposition,SVD) 是另一種正交矩陣分解法；SVD是最可靠的分解法，但是它比QR 分解法要花上近十倍的計算時間。[U,S,V]=svd(A)，其中U和V分別代表兩個正交矩陣，而S代表一對角矩陣。 和QR分解法相同， 原矩陣A不必為正方矩陣。使用SVD分解法的用途是解最小平方誤差法和數據壓縮。
 

矩陣的LL^T分解 

Cholesky 分解是把一個對稱正定的矩陣表示成一個下三角矩陣L和其轉置的乘積的分解。它要求矩陣的所有特征值必須大於零，故分解的下三角的對角元也是大於零的（LU三角分解法的變形）。

 

LDLT分解

對於對稱正定矩陣A，可以將A分解為A=LDL^T，並且該分解是唯一的，其中L為一下三角形單位矩陣（即主對角線元素皆為1），D為一對角矩陣（只在主對角線上有元素且不為0，其余皆為零），L^T為L的轉置矩陣。LDL^T 是Cholesky分解的變形：

LDL^T分解法實際上是Cholesky分解法的改進，因為Cholesky分解法雖然不需要選主元，但其運算過程中涉及到開方問題，而LDL^T分解法則避免了這一問題，可用於求解線性方程組。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    //線性方程求解 Ax =B;
    Matrix4d A;
    A << 2,-1,-1,1,
        1,1,-2,1,
        4,-6,2,-2,
        3,6,-9,7;

    Vector4d b(2,4,4,9);

    Vector4d x1 = A.colPivHouseholderQr().solve(b);
   

Vector4d x2

 = A.ldlt().solve(b));  // A sym. p.s.d.    #include <Eigen/Cholesky>

Vector4d x3

 = A.llt().solve(b));  // A sym. p.d.      #include <Eigen/Cholesky>

Vector4d x4

 = A.lu().solve(b));  // Stable and fast. #include <Eigen/LU>

Vector4d x5

 = A.qr().solve(b));  // No pivoting.     #include <Eigen/QR>

Vector4d x6

 = A.svd().solve(b));  // Stable, slowest. #include <Eigen/SVD>

// .ldlt() -> .matrixL() and .matrixD()
// .llt()  -> .matrixL()
// .lu()   -> .matrixL() and .matrixU()
// .qr()   -> .matrixQ() and .matrixR()
// .svd()  -> .matrixU(), .singularValues(), and .matrixV()

    std::cout << "The solution is:\n" << x <<"\n\n"<<x2<<"\n\n"<<x3 <<std::endl;
}

 

輸出為： 

 0
-1
-4
-3

-289.143
 448.714
 29.9082
 3.97959

  1.52903
   0.1758
-0.340206
0.0423223

 

結果不一致，這是為什么？

 

14、特征值分解

最小二乘法求解

最小二乘求解有兩種方式，jacobiSvd或者colPivHouseholderQr，4*4以下的小矩陣速度沒有區別，jacobiSvd可能更快，大矩陣最好用colPivHouseholderQr

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    MatrixXf A1 = MatrixXf::Random(3, 2);
    std::cout << "Here is the matrix A:\n" << A1 << std::endl;
    VectorXf b1 = VectorXf::Random(3);
    std::cout << "Here is the right hand side b:\n" << b1 << std::endl;
    //jacobiSvd 方式:Slow (but fast for small matrices)
    std::cout << "The least-squares solution is:\n"
    << A1.jacobiSvd(ComputeThinU | ComputeThinV).solve(b1) << std::endl;
    //colPivHouseholderQr方法:fast
    std::cout << "The least-squares solution is:\n"
    << A1.colPivHouseholderQr().solve(b1) << std::endl;
}

輸出為：

Here is the matrix A:
 0.680375   0.59688
-0.211234  0.823295
 0.566198 -0.604897
Here is the right hand side b:
-0.329554
 0.536459
-0.444451
The least-squares solution is:
-0.669626
 0.314253
The least-squares solution is:
-0.669626
 0.314253

 

15、稀疏矩陣

稀疏矩陣的頭文件包括：

#include <Eigen/SparseCore>
#include <Eigen/SparseCholesky>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/Sparse>

初始化有兩種方式：

1.使用三元組插入

typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;
triplets.reserve(estimation_of_entries); //estimation_of_entries是預估的條目
for(...)
{
    tripletList.push_back(T(i,j,v_ij));//第 i,j個有值的位置的值
}
SparseMatrixType mat(rows,cols);
mat.setFromTriplets(tripletList.begin(), tripletList.end());
// mat is ready to go!

2.直接將已知的非0值插入

SparseMatrix<double> mat(rows,cols);
mat.reserve(VectorXi::Constant(cols,6));
for(...)
{
    // i,j 個非零值 v_ij != 0
    mat.insert(i,j) = v_ij;
}
mat.makeCompressed(); // optional

稀疏矩陣支持大部分一元和二元運算:

sm1.real() sm1.imag() -sm1 0.5*sm1
sm1+sm2 sm1-sm2 sm1.cwiseProduct(sm2)

二元運算中，稀疏矩陣和普通矩陣可以混合使用

//dm表示普通矩陣
dm2 = sm1 + dm1;

16、C++數組和矩陣轉換

 使用Map函數，可以實現Eigen的矩陣和c++中的數組直接轉換,語法如下：

//@param MatrixType 矩陣類型
//@param MapOptions 可選參數，指的是指針是否對齊，Aligned, or Unaligned. The default is Unaligned.
//@param StrideType 可選參數，步長
/*
    Map<typename MatrixType,
        int MapOptions,
        typename StrideType>
*/
    int i;
    //數組轉矩陣
    double *aMat = new double[20];
    for(i =0;i<20;i++)
    {
        aMat[i] = rand()%11;
    }
    //靜態矩陣，編譯時確定維數 Matrix<double,4,5> 
    Eigen:Map<Matrix<double,4,5> > staMat(aMat);
 
 
    //輸出
    for (int i = 0; i < staMat.size(); i++)
        std::cout << *(staMat.data() + i) << " ";
    std::cout << std::endl << std::endl;
 
 
    //動態矩陣，運行時確定 MatrixXd
    Map<MatrixXd> dymMat(aMat,4,5);
 
 
    //輸出，應該和上面一致
    for (int i = 0; i < dymMat.size(); i++)
        std::cout << *(dymMat.data() + i) << " ";
    std::cout << std::endl << std::endl;
 
    //Matrix為列優先，如下返回指針
    dymMat.data();

 

 

 

矩陣的應用實例：

矩陣的定義 

#include<Eigen/Dense>

Matrix A; // Fixed rows and cols. Same as Matrix3d. Matrix B; // Fixed rows, dynamic cols.

Matrix C; // Full dynamic. Same as MatrixXd.

Matrix E; // Row major; default is column-major. Matrix3f P, Q, R; // 3x3 float matrix.

Vector3f x, y, z; // 3x1 float matrix.

RowVector3f a, b, c; // 1x3 float matrix.

VectorXd v; // Dynamic column vector of doubles

 Eigen的基礎使用

// Basic usage

// Eigen // Matlab // comments

x.size() // length(x) // vector size

C.rows() // size(C,1) // number of rows

C.cols() // size(C,2) // number of columns

x(i) // x(i+1) // Matlab is 1-based

C(i, j) // C(i+1,j+1) //

A.resize(4, 4); // Runtime error if assertions are on.

B.resize(4, 9); // Runtime error if assertions are on.

A.resize(3, 3); // Ok; size didn't change.

B.resize(3, 9); // Ok; only dynamic cols changed.

A <<1, 2, 3, // Initialize A. The elements can also be

4, 5, 6, // matrices, which are stacked along cols

7, 8, 9; // and then the rows are stacked.

B << A, A, A; // B is three horizontally stacked A's.

A.fill(10); // Fill A with all 10's.

 Eigen特殊矩陣生成

// Eigen // Matlab

MatrixXd::Identity(rows,cols) // eye(rows,cols)

C.setIdentity(rows,cols) // C = eye(rows,cols)

MatrixXd::Zero(rows,cols) // zeros(rows,cols)

C.setZero(rows,cols) // C = ones(rows,cols)

MatrixXd::Ones(rows,cols) // ones(rows,cols)

C.setOnes(rows,cols) // C = ones(rows,cols)

MatrixXd::Random(rows,cols) // rand(rows,cols)*2-1 // MatrixXd::Random returns uniform random numbers in (-1, 1).

C.setRandom(rows,cols) // C = rand(rows,cols)*2-1

VectorXd::LinSpaced(size,low,high) // linspace(low,high,size)'

v.setLinSpaced(size,low,high) // v = linspace(low,high,size)'

 Eigen矩陣分塊

// Matrix slicing and blocks. All expressions listed here are read/write.

// Templated size versions are faster. Note that Matlab is 1-based (a size N // vector is x(1)...x(N)).

// Eigen // Matlab

x.head(n) // x(1:n)

x.head() // x(1:n)

x.tail(n) // x(end - n + 1: end)

x.tail() // x(end - n + 1: end)

x.segment(i, n) // x(i+1 : i+n)

x.segment(i) // x(i+1 : i+n)

P.block(i, j, rows, cols) // P(i+1 : i+rows, j+1 : j+cols)

P.block(i, j) // P(i+1 : i+rows, j+1 : j+cols)

P.row(i) // P(i+1, :)

P.col(j) // P(:, j+1)

P.leftCols() // P(:, 1:cols)

P.leftCols(cols) // P(:, 1:cols)

P.middleCols(j) // P(:, j+1:j+cols)

P.middleCols(j, cols) // P(:, j+1:j+cols)

P.rightCols() // P(:, end-cols+1:end)

P.rightCols(cols) // P(:, end-cols+1:end)

P.topRows() // P(1:rows, :)

P.topRows(rows) // P(1:rows, :)

P.middleRows(i) // P(i+1:i+rows, :)

P.middleRows(i, rows) // P(i+1:i+rows, :)

P.bottomRows() // P(end-rows+1:end, :)

P.bottomRows(rows) // P(end-rows+1:end, :)

P.topLeftCorner(rows, cols) // P(1:rows, 1:cols)

P.topRightCorner(rows, cols) // P(1:rows, end-cols+1:end)

P.bottomLeftCorner(rows, cols) // P(end-rows+1:end, 1:cols)

P.bottomRightCorner(rows, cols) // P(end-rows+1:end, end-cols+1:end) P.topLeftCorner<rows,cols>() // P(1:rows, 1:cols)

P.topRightCorner<rows,cols>() // P(1:rows, end-cols+1:end)

P.bottomLeftCorner<rows,cols>() // P(end-rows+1:end, 1:cols)

P.bottomRightCorner<rows,cols>() // P(end-rows+1:end, end-cols+1:end)

 Eigen矩陣元素交換

// Of particular note is Eigen's swap function which is highly optimized. 

// Eigen 　　　　　　　　　　// Matlab

R.row(i) = P.col(j); 　　// R(i, :) = P(:, i)

R.col(j1).swap(mat1.col(j2)); // R(:, [j1 j2]) = R(:, [j2, j1])

Eigen矩陣轉置

// Views, transpose, etc; all read-write except for .adjoint().

// Eigen // Matlab

R.adjoint() // R'

R.transpose() // R.' or conj(R')

R.diagonal() // diag(R)

x.asDiagonal() // diag(x)

R.transpose().colwise().reverse(); // rot90(R)

R.conjugate() // conj(R)

 Eigen矩陣乘積

// All the same as Matlab, but matlab doesn't have *= style operators.

// Matrix-vector. Matrix-matrix.Matrix-scalar.

y = M*x; R = P*Q; R = P*s;

a = b*M; R = P - Q; R = s*P;

a *= M; R = P + Q; R = P/s;

R *= Q; R = s*P;

R += Q; R *= s;

R -= Q; R /= s;

 

Eigen矩陣單個元素操作

// Vectorized operations on each element independently // Eigen // Matlab

R = P.cwiseProduct(Q); // R = P .* Q

R = P.array() * s.array();// R = P .* s

R = P.cwiseQuotient(Q); // R = P ./ Q

R = P.array() / Q.array();// R = P ./ Q

R = P.array() + s.array();// R = P + s

R = P.array() - s.array();// R = P - s

R.array() += s; // R = R + s

R.array() -= s; // R = R - s

R.array()

R.array() <= Q.array(); // R <= Q

R.cwiseInverse(); // 1 ./ P

R.array().inverse(); // 1 ./ P

R.array().sin() // sin(P)

R.array().cos() // cos(P)

R.array().pow(s) // P .^ s

R.array().square() // P .^ 2

R.array().cube() // P .^ 3

R.cwiseSqrt() // sqrt(P)

R.array().sqrt() // sqrt(P)

R.array().exp() // exp(P)

R.array().log() // log(P)

R.cwiseMax(P) // max(R, P)

R.array().max(P.array()) // max(R, P)

R.cwiseMin(P) // min(R, P)

R.array().min(P.array()) // min(R, P)

R.cwiseAbs() // abs(P)

R.array().abs() // abs(P)

R.cwiseAbs2() // abs(P.^2)

R.array().abs2() // abs(P.^2)

(R.array() < s).select(P,Q); // (R < s ? P : Q)

 

Eigen矩陣簡化 

// Reductions.

int r, c;

// Eigen // Matlab

R.minCoeff() // min(R(:))

R.maxCoeff() // max(R(:))

s = R.minCoeff(&r, &c) // [s, i] = min(R(:)); [r, c] = ind2sub(size(R), i);

s = R.maxCoeff(&r, &c) // [s, i] = max(R(:)); [r, c] = ind2sub(size(R), i);

R.sum() // sum(R(:))

R.colwise().sum() // sum(R)

R.rowwise().sum() // sum(R, 2) or sum(R')'

R.prod() // prod(R(:))

R.colwise().prod() // prod(R)

R.rowwise().prod() // prod(R, 2) or prod(R')'

R.trace() // trace(R)

R.all() // all(R(:))

R.colwise().all() // all(R)

R.rowwise().all() // all(R, 2)

R.any() // any(R(:))

R.colwise().any() // any(R)

R.rowwise().any() // any(R, 2)

Eigen矩陣點乘

// Dot products, norms, etc.

// Eigen // Matlab
x.norm() // norm(x). Note that norm(R) doesn't work in Eigen. 
x.squaredNorm() // dot(x, x) Note the equivalence is not true for complex
x.dot(y) // dot(x, y)
x.cross(y) // cross(x, y) Requires #include

矩陣類型轉換

//// Type conversion

// Eigen // Matlab

A.cast(); // double(A)

A.cast(); // single(A)

A.cast(); // int32(A)

A.real(); // real(A)

A.imag(); // imag(A)

// if the original type equals destination type, no work is done

  

Eigen求解線性方程組Ax = b

// Solve Ax = b. Result stored in x. Matlab: x = A \ b.
x = A.ldlt().solve(b)); // A sym. p.s.d. #include 

x = A.llt() .solve(b)); // A sym. p.d. #include 

x = A.lu() .solve(b)); // Stable and fast. #include

x = A.qr() .solve(b)); // No pivoting. #include

x = A.svd() .solve(b)); // Stable, slowest. #include

// .ldlt() -> .matrixL() and .matrixD()

// .llt() -> .matrixL()

// .lu() -> .matrixL() and .matrixU()

// .qr() -> .matrixQ() and .matrixR()

// .svd() -> .matrixU(), .singularValues(), and .matrixV()

 Eigen矩陣特征值 

// Eigenvalue problems

// Eigen // Matlab

A.eigenvalues(); // eig(A);

EigenSolvereig(A); // [vecval] = eig(A)

eig.eigenvalues(); // diag(val)

eig.eigenvectors(); // vec

// For self-adjoint matrices use SelfAdjointEigenSolver<>