Eigen笔记之matrix

定义和初始化

matrix的原型为

1
2
3
4
5
6


Matrix<typename Scalar,
       int RowsAtCompileTime,
       int ColsAtCompileTime,
       int Options = 0,  // 代表存储方式是column-major(默认) 还是 row-major. 可使用常量 {RowMajor, ColMajor}
       int MaxRowsAtCompileTime = RowsAtCompileTime, //动态分配时行的上限
       int MaxColsAtCompileTime = ColsAtCompileTime>//动态分配时列的上限

我们可以定义fixed和dynamic的matrix如下

1
2
3
4
5


typedef Matrix<float, 4, 4> Matrix4f;
typedef Matrix<float, 3, 1> Vector3f;
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi;
typdef Matrix<float, 3, Dynamic> Matrix3Xf;

声明：

1
2
3
4
5
6
7
8
9


// a is a 3-by-3 matrix, with a plain float[9] array of
// uninitialized coefficients,
Matrix3f a;
// b is a dynamic-size matrix whose size is currently 0-by-0,
// and whose array of coefficients hasn't yet been allocated at all.
MatrixXf b;

MatrixXf a(10,15);
VectorXf b(30);

可以使用,的方式进行初始化：

1
2
3
4
5


Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;
std::cout << m;

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


MatrixXf matA(2, 2);
matA << 1, 2, 3, 4;
MatrixXf matB(4, 4);
matB << matA, matA/10, matA/10, matA;
std::cout << matB << std::endl;
=>
1   2   0.1 0.2
3   4   0.3 0.4
0.1 0.2   1   2
0.3 0.4   3   4

1
2
3
4
5
6
7
8
9


Matrix3f m;
m.row(0) << 1, 2, 3;
m.block(1,0,2,2) << 4, 5, 7, 8;
m.col(2).tail(2) << 6, 9;
std::cout << m;
=>
1 2 3
4 5 6
7 8 9

基本操作

元素获取

通过使用(i, j)进行元素的获取

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18


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

如果使用（i）将是访问第0列的第i位元素 (取决于matrx的存储方式是column-major还是row-major，默认为column-major)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


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);
cout << m(0) << endl;
cout << m(1) << endl;
=>
3
2.5

获取某行或者某列使用matrix.row(i),matrix.col(i):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


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

shape相关

rows() 获取行数
cols() 获取列数
size() 获取元素的个数
resize() 修改shape (元素的值可能会产生变化)
conservativeResize() 修改shape (元素的值可保留)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

int main()
{
    MatrixXd m(2,5);
    m.resize(4,3);
    std::cout << "The matrix m is of size "
          << m.rows() << "x" << m.cols() << std::endl;
    std::cout << "It has " << m.size() << " coefficients" << std::endl;
    VectorXd v(2);
    v.resize(5);
    std::cout << "The vector v is of size " << v.size() << std::endl;
    std::cout << "As a matrix, v is of size "

赋值

setZero()
fill(val)
setOnes()
setRandom()

1
2
3
4
5
6
7
8
9


  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;
  m1.setZero(); // 赋值全为0
  m1.fill(2);  // 赋值全为2
  m1.setOnes(); // 赋值全为1
  m1.setRandom(); // 赋值为随机数

获取数组指针

使用函数data()即可。在底层matrix和vector都是基于一维数组实现的,所以调用该函数之后得到的即是一个一维数组指针Type* prt(需要注意默认存储方式是column-major).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29


#include <iostream>
#include <eigen/Eigen/Dense>

using Eigen::Dynamic;
typedef Eigen::Matrix<double, Dynamic, Dynamic> MatrixD;

using namespace std;

int main() {
    MatrixD X2(4, 4);
    double cnt = 1;
    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++) {
            X2(i, j) = cnt++;
        }
    }
    double* arr = X2.data();
    cout << X2 << endl;
    for (int i = 0; i < 16; i++) {
        cout << arr[i] << " ";
    }
    return 0;
}
=>
1  2  3  4
5  6  7  8
9 10 11 12
13 14 15 16
1 5 9 13 2 6 10 14 3 7 11 15 4 8 12 16

算术运算

加减运算
- binary operator + as in a+b
- binary operator - as in a-b
- unary operator - as in -a
- compound operator += as in a+=b
- compound operator -= as in a-=b
与标量的乘除运算
- binary operator * as in matrix*scalar
- binary operator * as in scalar*matrix
- binary operator / as in matrix/scalar
- compound operator *= as in matrix *=scalar
- compound operator /= as in matrix/=scalar

需要注意在Eigen库中没有实现matrix与标量进行加减的运算！所以我们需要自己写个函数来实现加减运算，比如：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14


  // Eigen Matrix plus a scalar
  template<class Derived, class Type>
  void PlusScalar (MatrixBase<Derived>& m1, Type num) {
    int rows = m1.rows(), cols = m1.cols();
    m1 = m1 + num * MatrixBase<Derived>::Constant(rows, cols, num);
  }

  // Eigen Matrix minus a scalar
  template<class Derived, class Type>
  void MinusScalar (MatrixBase<Derived>& m1, Type num) {
    int rows = m1.rows(), cols = m1.cols();
    m1 = m1 - num * MatrixBase<Derived>::Constant(rows, cols, num);
  }

转置与共轭

transpose()
conjugate()
adjoint()

1
2
3
4
5


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;

In place 的转置

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


Matrix2i a; a << 1, 2, 3, 4;
cout << "Here is the matrix a:\n" << a << endl;
a = a.transpose(); // !!! do NOT do this !!!
cout << "and the result of the aliasing effect:\n" << a << endl;
=>
Here is the matrix a:
1 2
3 4
and the result of the aliasing effect:
1 2
2 4

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;

矩阵与矩阵，矩阵与向量之间的乘法

binary operator * as in a*b
compound operator *= as in a *=b 执行点乘运算

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37


#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

注:为了避免aliasing的问题，可以使用noalias()

1

c.noalias() += a * b;

点乘(dot produc)和叉乘运算(cross product)

dot()
corss()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22


#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

Reduction operations

sum()
prod()
maxCoeff()
minCoeff()
trace() 计算对角线元素的和

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23


#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  cout << "Here is mat.sum():       " << mat.sum()       << endl;
  cout << "Here is mat.prod():      " << mat.prod()      << endl;
  cout << "Here is mat.mean():      " << mat.mean()      << endl;
  cout << "Here is mat.minCoeff():  " << mat.minCoeff()  << endl;
  cout << "Here is mat.maxCoeff():  " << mat.maxCoeff()  << endl;
  cout << "Here is mat.trace():     " << mat.trace()     << endl;
}
=>
Here is mat.sum():       10
Here is mat.prod():      24
Here is mat.mean():      2.5
Here is mat.minCoeff():  1
Here is mat.maxCoeff():  4
Here is mat.trace():     5

Norm computations

squaredNorm() 求向量v的$||v||^2$
norm() 求向量的$||V||$
lpNorm<p>() p-norm

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37


#include <Eigen/Dense>
#include <iostream>

using namespace std;
using namespace Eigen;

int main()
{
  VectorXf v(2);
  MatrixXf m(2,2), n(2,2);

  v << -1,
       2;

  m << 1,-2,
       -3,4;

  cout << "v.squaredNorm() = " << v.squaredNorm() << endl;
  cout << "v.norm() = " << v.norm() << endl;
  cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << endl;
  cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << endl;

  cout << endl;
  cout << "m.squaredNorm() = " << m.squaredNorm() << endl;
  cout << "m.norm() = " << m.norm() << endl;
  cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << endl;
}
=>
v.squaredNorm() = 5
v.norm() = 2.23607
v.lpNorm<1>() = 3
v.lpNorm<Infinity>() = 2

m.squaredNorm() = 30
m.norm() = 5.47723
m.lpNorm<1>() = 10
m.lpNorm<Infinity>() = 4

Boolean reductions

all() 检查是否所有的元素都为true
any() 检查是否存在至少一个元素为true
count() 计算元素中true的个数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19


#include <Eigen/Dense>
#include <iostream>

using namespace std;
using namespace Eigen;

int main()
{
  ArrayXXf a(2,2);
  a << 1,2,
       3,4;
  cout << "(a > 0).all()   = " << (a > 0).all() << endl;
  cout << "(a > 0).any()   = " << (a > 0).any() << endl;
  cout << "(a > 0).count() = " << (a > 0).count() << endl;
  cout << endl;
  cout << "(a > 2).all()   = " << (a > 2).all() << endl;
  cout << "(a > 2).any()   = " << (a > 2).any() << endl;
  cout << "(a > 2).count() = " << (a > 2).count() << endl;
}

Visitors

maxCoeff(&x, &y)
minCoeff(&x, &y)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29


#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
  Eigen::MatrixXf m(2,2);

  m << 1, 2,
       3, 4;

  //get location of maximum
  MatrixXf::Index maxRow, maxCol;
  float max = m.maxCoeff(&maxRow, &maxCol);

  //get location of minimum
  MatrixXf::Index minRow, minCol;
  float min = m.minCoeff(&minRow, &minCol);

  cout << "Max: " << max <<  ", at: " <<
     maxRow << "," << maxCol << endl;
  cout << "Min: " << min << ", at: " <<
     minRow << "," << minCol << endl;
}
=>
Max: 4, at: 1,1
Min: 1, at: 0,0

Partial reductions

colwise()
rowwise()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;

  std::cout << "Column's maximum: " << std::endl
   << mat.colwise().maxCoeff() << std::endl;
}
=>
Column's maximum:
3 2 7 9

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;

  std::cout << "Row's maximum: " << std::endl
   << mat.rowwise().maxCoeff() << std::endl;
}
=>
Row's maximum:
9
7

特殊矩阵和向量

Zero()
Constant(value)
Random() 对于浮点数生成的随机数范围是[-1, 1]，对于正数生成的随机数范围是整个取值空间
Identity()
LinSpaced(size, low, high)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


std::cout << "A fixed-size array:\n";
Array33f a1 = Array33f::Zero();
std::cout << a1 << "\n\n";

std::cout << "A one-dimensional dynamic-size array:\n";
ArrayXf a2 = ArrayXf::Zero(3);
std::cout << a2 << "\n\n";

std::cout << "A two-dimensional dynamic-size array:\n";
ArrayXXf a3 = ArrayXXf::Zero(3, 4);
std::cout << a3 << "\n";

MatrixD ones = MatrixD::Constant(2, 2, 0);

基于向量生成对角矩阵asDiagonal()

1
2
3
4
5
6


  typedef Eigen::Matrix<double, 1, Dynamic> VectorD;
  VectorD V3(4);
  V3 << 2, 3, 4, 5;
  MatrixD M2;
  M2 = V3.asDiagonal() * M;

Map

将c/c++中的array转化为Eigen中的vector和matrix。

Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >

初始化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


Map<MatrixXf> mf(pf,rows,columns);
Map<const Vector4i> mi(pi);

// convert STL vector to Eigen Vector
vector<int> local_vec;
for (int i = 0; i < 3; i++) {
    local_vec.push_back(i);
}
Map<MatrixXi> local_v(local_vec.data(), 1, local_vec.size());

vector<double> mv;
for (int i = 0; i < 5; i++) {
  for (int j = 0; j < 3; j++) {
    mv.push_back(i);
  }
}
Map<MatrixD> X(mv.data(), 3, 5);

使用

Map对象可以像Eigen中的其它类型一样使用。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28


typedef Matrix<float,1,Dynamic> MatrixType;
typedef Map<MatrixType> MapType;
typedef Map<const MatrixType> MapTypeConst;   // a read-only map
const int n_dims = 5;

MatrixType m1(n_dims), m2(n_dims);
m1.setRandom();
m2.setRandom();
float *p = &m2(0);  // get the address storing the data for m2
MapType m2map(p,m2.size());   // m2map shares data with m2
MapTypeConst m2mapconst(p,m2.size());  // a read-only accessor for m2

cout << "m1: " << m1 << endl;
cout << "m2: " << m2 << endl;
cout << "Squared euclidean distance: " << (m1-m2).squaredNorm() << endl;
cout << "Squared euclidean distance, using map: " <<
  (m1-m2map).squaredNorm() << endl;
m2map(3) = 7;   // this will change m2, since they share the same array
cout << "Updated m2: " << m2 << endl;
cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;
/* m2mapconst(2) = 5; */   // this yields a compile-time error
=>
m1:   0.68 -0.211  0.566  0.597  0.823
m2: -0.605  -0.33  0.536 -0.444  0.108
Squared euclidean distance: 3.26
Squared euclidean distance, using map: 3.26
Updated m2: -0.605  -0.33  0.536      7  0.108
m2 coefficient 2, constant accessor: 0.536

改变Map的数组需要使用new:

1
2
3
4
5
6
7
8


int data[] = {1,2,3,4,5,6,7,8,9};
Map<RowVectorXi> v(data,4);
cout << "The mapped vector v is: " << v << "\n";
new (&v) Map<RowVectorXi>(data+4,5);
cout << "Now v is: " << v << "\n";
=>
The mapped vector v is: 1 2 3 4
Now v is: 5 6 7 8 9

1
2
3
4
5
6
7


Map<Matrix3f> A(NULL);  // don't try to use this matrix yet!
VectorXf b(n_matrices);
for (int i = 0; i < n_matrices; i++)
{
  new (&A) Map<Matrix3f>(get_matrix_pointer(i));
  b(i) = A.trace();
}

Reshape

Eigen中没有方便的方法用来执行切分和reshape一个matrix,要执行这个功能需要借助Map类。

reshape

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


MatrixXf M1(3,3);    // Column-major storage
M1 << 1, 2, 3,
      4, 5, 6,
      7, 8, 9;

Map<RowVectorXf> v1(M1.data(), M1.size());
cout << "v1:" << endl << v1 << endl;

Matrix<float,Dynamic,Dynamic,RowMajor> M2(M1);
Map<RowVectorXf> v2(M2.data(), M2.size());
cout << "v2:" << endl << v2 << endl;
=>
v1:
1 4 7 2 5 8 3 6 9
v2:
1 2 3 4 5 6 7 8 9

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


MatrixXf M1(2,6);    // Column-major storage
M1 << 1, 2, 3,  4,  5,  6,
      7, 8, 9, 10, 11, 12;

Map<MatrixXf> M2(M1.data(), 6,2);
cout << "M2:" << endl << M2 << endl;
=>
M2:
 1  4
 7 10
 2  5
 8 11
 3  6
 9 12

Slicing

1
2
3
4
5
6
7
8


RowVectorXf v = RowVectorXf::LinSpaced(20,0,19);
cout << "Input:" << endl << v << endl;
Map<RowVectorXf,0,InnerStride<2> > v2(v.data(), v.size()/2);
cout << "Even:" << v2 << endl;
=>
Input:
 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
Even: 0  2  4  6  8 10 12 14 16 18

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28


MatrixXf M1 = MatrixXf::Random(3,8);
cout << "Column major input:" << endl << M1 << "\n";
Map<MatrixXf,0,OuterStride<> > M2(M1.data(), M1.rows(), (M1.cols()+2)/3, OuterStride<>(M1.outerStride()*3));
cout << "1 column over 3:" << endl << M2 << "\n";

typedef Matrix<float,Dynamic,Dynamic,RowMajor> RowMajorMatrixXf;
RowMajorMatrixXf M3(M1);
cout << "Row major input:" << endl << M3 << "\n";
Map<RowMajorMatrixXf,0,Stride<Dynamic,3> > M4(M3.data(), M3.rows(), (M3.cols()+2)/3,
Stride<Dynamic,3>(M3.outerStride(),3));
cout << "1 column over 3:" << endl << M4 << "\n";
=>
Column major input:
   0.68   0.597   -0.33   0.108   -0.27   0.832  -0.717  -0.514
 -0.211   0.823   0.536 -0.0452  0.0268   0.271   0.214  -0.726
  0.566  -0.605  -0.444   0.258   0.904   0.435  -0.967   0.608
1 column over 3:
   0.68   0.108  -0.717
 -0.211 -0.0452   0.214
  0.566   0.258  -0.967
Row major input:
   0.68   0.597   -0.33   0.108   -0.27   0.832  -0.717  -0.514
 -0.211   0.823   0.536 -0.0452  0.0268   0.271   0.214  -0.726
  0.566  -0.605  -0.444   0.258   0.904   0.435  -0.967   0.608
1 column over 3:
   0.68   0.108  -0.717
 -0.211 -0.0452   0.214
  0.566   0.258  -0.967

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25


  MatrixD m1 = MatrixD::Random(8, 8);
  int cnt = 0;
  for (int i = 0; i < 8; i++) {
    for (int j = 0; j < 8; j++) {
      m1(i, j) = cnt++;
    }
  }
  cout << m1 << endl;
  Map<MatrixD, 0, Stride<Dynamic, Dynamic> > m2(m1.data(), 3, 3, Stride<Dynamic, Dynamic>(m1.rows(), 2));
  // 其中的Dynamic也可以是具体的数值
  cout << m2 << endl;
==>
0  1  2  3  4  5  6  7
8  9 10 11 12 13 14 15
16 17 18 19 20 21 22 23
24 25 26 27 28 29 30 31
32 33 34 35 36 37 38 39
40 41 42 43 44 45 46 47
48 49 50 51 52 53 54 55
56 57 58 59 60 61 62 63


0  2  4
8 10 12
16 18 20

其中Stride<>的定义为Stride<OuterStride, InnerStride>, OuterStide代表组间距离步长(一般为行数或者列数)，InnerStride代表组内步长

Map的原型为Map(dataPtr, rows, cols, StrideType())，rows和cols代表的生成的矩阵大小

In place function

MatrixBase::adjointInPlace()
DenseBase::reverseInPlace()
DenseBase::transposeInPlace()

Aliasing

Aliasing问题会出现在同一个matrix或者array出现在赋值语句的两边，比如m=m.transpose()

可以使用eval()函数来避免该问题

1
2
3
4
5
6
7


MatrixXi mat(3,3);
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;

// The eval() solves the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();
cout << "After the assignment, mat = \n" << mat << endl;

Storage orders

矩阵在内存中存储有两种方式row-major和column-major.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


Matrix<int, 3, 4, ColMajor> Acolmajor;
Acolmajor << 8, 2, 2, 9,
             9, 1, 4, 4,
             3, 5, 4, 5;
cout << "The matrix A:" << endl;
cout << Acolmajor << endl << endl;

cout << "In memory (column-major):" << endl;
for (int i = 0; i < Acolmajor.size(); i++)
  cout << *(Acolmajor.data() + i) << "  ";
cout << endl << endl;

Matrix<int, 3, 4, RowMajor> Arowmajor = Acolmajor;
cout << "In memory (row-major):" << endl;
for (int i = 0; i < Arowmajor.size(); i++)
  cout << *(Arowmajor.data() + i) << "  ";
cout << endl;

将Eigen类型作为函数的参数

使用模板

EigenBase 示例：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18


#include <iostream>
#include <Eigen/Core>
using namespace Eigen;

template <typename Derived>
void print_size(const EigenBase<Derived>& b)
{
  std::cout << "size (rows, cols): " << b.size() << " (" << b.rows()
            << ", " << b.cols() << ")" << std::endl;
}

int main()
{
    Vector3f v;
    print_size(v);
    // v.asDiagonal() returns a 3x3 diagonal matrix pseudo-expression
    print_size(v.asDiagonal());
}

DenseBase示例：

1
2
3
4
5


template <typename Derived>
void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
{
  std::cout << "block: " << b.block(x,y,r,c) << std::endl;
}

ArrayBase示例：

1
2
3
4
5


template <typename Derived>
void print_max_coeff(const ArrayBase<Derived> &a)
{
  std::cout << "max: " << a.maxCoeff() << std::endl;
}

MatrixBase示例：

1
2
3
4
5
6
7


template <typename Derived>
void print_inv_cond(const MatrixBase<Derived>& a)
{
  const typename JacobiSVD<typename Derived::PlainObject>::SingularValuesType&
  sing_vals = a.jacobiSvd().singularValues();
  std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
}

Multiple templated arguments 示例：

1
2
3
4
5


template <typename DerivedA,typename DerivedB>
typename DerivedA::Scalar squaredist(const MatrixBase<DerivedA>& p1,const MatrixBase<DerivedB>& p2)
{
  return (p1-p2).squaredNorm();
}

Block operations

block operations提供对matrix进行整块数据的操作。

1

matrix.block(i,j,p,q);  // Block of size (p,q), starting at (i,j)

Block operation	Version constructing a dynamic-size block expression	Version constructing a fixed-size block expression
Top-left p by q block *	`matrix.topLeftCorner(p,q);`	`matrix.topLeftCorner<p,q>();`
Bottom-left p by q block *	`matrix.bottomLeftCorner(p,q);`	`matrix.bottomLeftCorner<p,q>();`
Top-right p by q block *	`matrix.topRightCorner(p,q);`	`matrix.topRightCorner<p,q>();`
Bottom-right p by q block *	`matrix.bottomRightCorner(p,q);`	`matrix.bottomRightCorner<p,q>();`
Block containing the first q rows *	`matrix.topRows(q);`	`matrix.topRows<q>();`
Block containing the last q rows *	`matrix.bottomRows(q);`	`matrix.bottomRows<q>();`
Block containing the first p columns *	`matrix.leftCols(p);`	`matrix.leftCols<p>();`
Block containing the last q columns *	`matrix.rightCols(q);`	`matrix.rightCols<q>();`

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33


#include <Eigen/Dense>
#include <iostream>

using namespace std;

int main()
{
  Eigen::Matrix4f m;
  m << 1, 2, 3, 4,
       5, 6, 7, 8,
       9, 10,11,12,
       13,14,15,16;
  cout << "m.leftCols(2) =" << endl << m.leftCols(2) << endl << endl;
  cout << "m.bottomRows<2>() =" << endl << m.bottomRows<2>() << endl << endl;
  m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();
  cout << "After assignment, m = " << endl << m << endl;
}
=>
m.leftCols(2) =
 1  2
 5  6
 9 10
13 14

m.bottomRows<2>() =
 9 10 11 12
13 14 15 16

After assignment, m =
 8 12 16  4
 5  6  7  8
 9 10 11 12
13 14 15 16

block operations提供对vector进行整块数据的操作。

Block operation	Version constructing a dynamic-size block expression	Version constructing a fixed-size block expression
Block containing the first n elements *	vector.head(n);	`vector.head<n>();`
Block containing the last n elements *	vector.tail(n);	`vector.tail<n>();`
Block containing n elements, starting at position i *	vector.segment(i,n);	`vector.segment<n>(i);`

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32


#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

常用函数

寻找矩阵中最小的值且得到该值的位置

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22


#include <iostream>
#include <eigen/Eigen/Dense>
#include <vector>

typedef Eigen::Matrix<double, Dynamic, Dynamic> MatrixD;
using namespace std;

int main() {
  MatrixD m1 (4, 4);
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 4; j++) {
      m1(i, j) = i + j + 8;
    }
  }
  m1(1, 2) = 4;
  cout << m1 << endl;
  int index1;
  int index2;
  m1.minCoeff(&index1, &index2);
  cout << index1 << " " << index2 << endl;
  return 0;
}

寻找矩阵中每一行的最小的值且得到该值的位置

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24


#include <iostream>
#include <eigen/Eigen/Dense>
#include <vector>

typedef Eigen::Matrix<double, Dynamic, Dynamic> MatrixD;
using namespace std;

int main() {
  MatrixD m1 (4, 4);
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 4; j++) {
      m1(i, j) = i + j + 8;
    }
  }
  m1(1, 2) = 4;
  cout << m1 << endl;
  int index1;
  for (int i = 0; i < 4; i++) {
    m1.row(i).minCoeff(&index1);
    cout << index1 << endl;
  }
  return 0;
}

对矩阵或向量中的所有元素求平方和

1
2
3
4
5
6


  MatrixD A(4, 4);
  A << 1, 2, 0, 0, 0, 1, 0, 0, 0, 0, 1, 2, 0, 0, 0, 1;
  double s = A.squaredNorm();
  cout << s << endl;
  cout << A << endl;   // 原矩阵保持不变

求是对称矩阵的最大k个特征值对应的特征向量

1
2
3
4
5


#include <Eigen/Dense>
...
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXf> solver(n);
solver.compute(L);
Eigen::MatrixXf X = solver.eigenvectors().rowwise().reverse().block(0,0,n,k);

高级模块

LU 分解

通过调用Eigen::FullPivLU<>之后会得到一个矩阵，该矩阵将L矩阵和U矩阵共同存储在一个矩阵之中(具体请看下面的例子)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47


typedef Matrix<double, 5, 3> Matrix5x3;
typedef Matrix<double, 5, 5> Matrix5x5;
Matrix5x3 m = Matrix5x3::Random();
cout << "Here is the matrix m:" << endl << m << endl;
Eigen::FullPivLU<Matrix5x3> lu(m);
cout << "Here is, up to permutations, its LU decomposition matrix:"
     << endl << lu.matrixLU() << endl;
cout << "Here is the L part:" << endl;
Matrix5x5 l = Matrix5x5::Identity();
l.block<5,3>(0,0).triangularView<StrictlyLower>() = lu.matrixLU();
cout << l << endl;
cout << "Here is the U part:" << endl;
Matrix5x3 u = lu.matrixLU().triangularView<Upper>();
cout << u << endl;
cout << "Let us now reconstruct the original matrix m:" << endl;
cout << lu.permutationP().inverse() * l * u * lu.permutationQ().inverse() << endl;
=>
Here is the matrix m:
   0.68  -0.605 -0.0452
 -0.211   -0.33   0.258
  0.566   0.536   -0.27
  0.597  -0.444  0.0268
  0.823   0.108   0.904
Here is, up to permutations, its LU decomposition matrix:
 0.904  0.823  0.108
-0.299  0.812  0.569
 -0.05  0.888   -1.1
0.0296  0.705  0.768
 0.285 -0.549 0.0436
Here is the L part:
     1      0      0      0      0
-0.299      1      0      0      0
 -0.05  0.888      1      0      0
0.0296  0.705  0.768      1      0
 0.285 -0.549 0.0436      0      1
Here is the U part:
0.904 0.823 0.108
    0 0.812 0.569
    0     0  -1.1
    0     0     0
    0     0     0
Let us now reconstruct the original matrix m:
   0.68  -0.605 -0.0452
 -0.211   -0.33   0.258
  0.566   0.536   -0.27
  0.597  -0.444  0.0268
  0.823   0.108   0.904

Cholesky 分解

标准的$LL^T$分解

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21


#include <iostream>
#include <Library/eigen/Eigen/Dense>

using namespace Eigen;

using namespace std;
int main(){
  MatrixXd A(3,3);
  A << 4,-1,2, -1,6,0, 2,0,5;
  cout << "The matrix A is" << endl << A << endl;

  LLT<MatrixXd> lltOfA(A); // compute the Cholesky decomposition of A
  MatrixXd L = lltOfA.matrixL(); // retrieve factor L  in the decomposition
  // The previous two lines can also be written as "L = A.llt().matrixL()"

  cout << "The Cholesky factor L is" << endl << L << endl;
  cout << "To check this, let us compute L * L.transpose()" << endl;
  cout << L * L.transpose() << endl;
  cout << "This should equal the matrix A" << endl;

}

利用Cholesky分解求解方程

求解方程$Ax=b$,使用该方法默认矩阵$A$可逆。

例子1：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


#include <iostream>
#include <Library/eigen/Eigen/Dense>

using namespace Eigen;
typedef Matrix<float,Dynamic,2> DataMatrix;
using namespace std;
int main(){
  // let's generate some samples on the 3D plane of equation z = 2x+3y (with some noise)
  DataMatrix samples = DataMatrix::Random(12,2);
  VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
  // and let's solve samples * [x y]^T = elevations in least square sense:
  Matrix<float,2,1> xy
    = (samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations));
  cout << xy << endl;
}

例子2：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23


#include <iostream>
#include <Library/eigen/Eigen/Dense>

using Eigen::Dynamic;
typedef Eigen::Matrix<double, Dynamic, Dynamic> MatrixD;
using namespace Eigen;

using namespace std;
int main(){
  MatrixD A(4, 4);
  A << 1, 2, 3, 4, 0, 1, 2, 3, 0, 0, 1, 2, 0, 0, 0, 1;
  MatrixD x(4, 1);
  MatrixD b(4, 1);
  b << 1, 2, 3, 4;
  // LLT<MatrixD> llt;
  // llt.compute(A.adjoint() * A);
  // x = llt.solve(A.adjoint() * b);
  x = (A.adjoint() * A).llt().solve(A.adjoint() * b);
  cout << A << endl;
  cout << "-------------" << endl;
  cout << b << endl;
  cout << "---------------" << endl;
  cout << x;

对角矩阵

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15


auto mat = vec.asDiagonal();
auto mat = vec.asDiagonal().toDenseMatrix();

DiagonalMatrix<Scalar,SizeAtCompileTime> diag1(size);
diag1.diagonal() = vector;

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
mat3  = scalar * diag1 * mat1;
mat3 += scalar * mat1 * vec1.asDiagonal();
mat3 = vec1.asDiagonal().inverse() * mat1
mat3 = mat1 * diag1.inverse()

特殊问题

将fixed的vector或matrix作为类的内部成员变量需要使用宏 EIGEN_MAKE_ALIGNED_OPERATOR_NEW

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


class Foo
{
    ...
    Eigen::Vector2d v;
    ...
    public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

 Foo *foo = new Foo;

访问元素必须是Int类型

1
2
3
4
5


  VectorD v(4);
  v << 1, 2, 3, 4;
  double n = 2;
  cout << v((int)n) << endl;
  cout << v(n) << endl;  // error

构造$[x_1^T, x_1^T, \cdots, x_1^T]$的矩阵

MatrixD all_one = MatrixD::Constant(out_features, 1, 1);
m = (all_one * x1).transpose();

$(\mathbf{1} * x_1)^T$其中$\mathbf{1}$是一个全为1的列向量

直接return 一个matrix 不会重新分配保存数据的内存

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


MatrixD test() {
  MatrixD a(2, 2);
  a << 1, 2, 3, 4;
  cout << "a address: " << a.data() << endl;
  return a;
}

int main() {
  MatrixD b =  test();
  cout << "b address: " << b.data() << endl;
  return 0;
}
// 输出的两个地址一样

cwise product of matrix

1
2
3
4
5
6


  MatrixD b(2, 2);
  MatrixD a(2, 2);
  b << 1, 2, 3, 4;
  a << 2, 3, 3, 3;
  auto c = a.cwiseProduct(b);
  auto c = a.cwiseQuotient(b);

Cast a matrix to different type

1
2


Eigen::MatrixXd d;                       // Matrix of doubles.
Eigen::MatrixXf f = d.cast <float> ();   // Matrix of floats.

修改后的Eigen，添加的函数resize(DataType*, int , int), reassign(DataType*).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28


#include "core/base.h"

using namespace std;

void func(MatrixD& x) {
  MatrixD y(4, 4);
  y.setOnes();
  MatrixD z = x * y;
  cout << "z:\n" << z << endl;
}
void Test() {
 double* ptr = (double*)malloc(16 * sizeof(double));
 for (int i = 0; i < 16; i++)
   ptr[i] = i;
 MatrixD x;
 x.resize(ptr, 4, 4);
 // MatrixD x(4, 4);
 // x.reassign(ptr);
 cout << "assigned: " << x << endl;
 x.setOnes();
 func(x);
}

int main(int argc, char ** argv) {
  Test();
  return 0;
}