eigen 筆記1

c++ 的 eigen 類似於 python 的 numpy, 還有一個類似的庫是 Armadillo, 當然還有 opencv.

Armadillo 與 matlab 在函數名稱上更接近, 但是 TensorFlow 和 Ceres 使用了 eigen.

這里不講究誰優誰劣, 入門階段迅速掌握一個, 用起來就夠了.

 

1. The Matrix Class

1) The first three template parameters of Matrix

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

typedef Matrix<float, 4, 4> Matrix4f;

 

2) Vectors

In Eigen, vectors are just a special case of matrices.

typedef Matrix<float, 3, 1> Vector3f;

typedef Matrix<int, 1, 2> RowVector2i;

 

3) The special value dynamic

Matrices dimensions can be unknown at compile time in Eigen.

typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

typedef Matrix<int, Dynamic, 1> VectorXi;

Matrix<float, 3, Dynamic>;

 

4) Constructors

Matrix3f a;

MatrixXf b;

MatrixXf a(10, 15);

MatrixXf b(30);

Matrix3f a(3, 3);

Vector2d a(5.0, 6.0);

Vector3d b(5.0, 6.0, 7.0);

Matrix3f 聲明的矩陣, 大小是固定的, 不能修改. 4維之下的矩陣可以使用形如 Matrix4f 的形式, 稍微增加編譯時間, 但是執行更快.

 

5) Coefficient accessors

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

m(index) 除了可以訪問 vector, 還可以訪問 matrix, matrix 在 Eigen 中默認按列存儲. 在上例中, m(2) == -1.

[] 也被重載了, 類似 () 的功能, 但是只能有一個 index, 因為 c++ 的特性, matrix[i, j] == matrix[j].

推測, 訪問單獨或批量的列或行, 應該有專有的函數, 用 () 和 [] 只能訪問單個值.

 

6) Comma-initialization

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

 

7) Resizing

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 "
            << v.rows() << "x" << v.cols() << std::endl;
}

resize() 會覆蓋之前矩陣的值, 如果要保留原有未知的值, 可以使用 conservativeResize().

再次強調, 形如 Matrix4d 聲明的矩陣大小不可更改.

 

8) Assignment and resizing

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;

如果等號左側是 fixed size, 則 resizing 不被允許.

 

9) Convenience typedefs

MatrixNt for Matrix<type, N, N>

VectorNt for Matrix<type, N, 1>

RowVectorNt for Matrix<type, 1, N>

 

N can be 2, 3,... or X (Dynamic)

t can be i (int), f (float), d (double), cf (complex<float>), cd (complex<double>)

 

2. Matrix and vector arithmetic

1) Introduction

Eigen 重載了 matrix 之間的運算符號, 但是 matrix 和 scalar 之間的運算符號未重載. 

 

2) Addition and subtraction

Eigen 不會對 scalar 進行自動類型轉換.

a + b, a - b,  -a, a += b, a -= b

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

 

3) Scalar multiplication and division

matrix * scalar, scalar * matrix, matrix / scalar, matrix *= scalar, matrix /= scalar

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

 

4) A note about expression templates

簡單說就是 Eigen 對 arithmetic expression 有進行優化, 但這是它應該做的.

 

5) Transposition and conjugation

transpose() 轉置

MatrixXcf a = MatrixXd::Random(2,2);
cout << "Here is the matrix a\n" << a << endl;
cout << "Here is the matrix a^T\n" << a.transpose() << endl;

不能使用 a = a.transpose(), a 的值會發生錯亂, 要使用 a = a.transposeInPlace().

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;
a.transposeInPlace();
cout << "and after being transposed:\n" << a << endl;

 

6) Matrix-matrix and matrix-vector multiplication

a * b, a *= b

對於 a *= b, 如果 size 會發生變化, 則 a 不能是 fixed size

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

 

7) Dot product and cross product

dot() 點積

cross() 叉積

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

 

8) Basic arithmetic reduction operations

sum(), prod(), mean(), minCoeff(), maxCoeff(), trace()

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

a.trace() == a.diagonal().sum()

minCoeff() 和 maxCoeff() 除了返回最小最大值, 還可以獲取該值的索引

  Matrix3f m = Matrix3f::Random();
  std::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 << ")\n\n";
  RowVector4i v = RowVector4i::Random();
  int maxOfV = v.maxCoeff(&i);
  cout << "Here is the vector v: " << v << endl;
  cout << "Its maximum coefficient (" << maxOfV 
       << ") is at position " << i << endl;

 

9) Validity of operations

編譯時和運行時都會檢查運算的合法性

Matrix3f m;
Vector4f v;
v = m*v;      // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES

MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"

 

3. The array class and coefficient-wise operations

1) Array types

Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

可以有多維數組

Array<float, Dynamic, 1> ArrayXf

Array<float, 3, 1> Array3f

Array<double, Dynamic, Dynamic> ArrayXXd

Array<double, 3, 3> Array33d

 

2) Accessing values inside an Array

int main()
{
  ArrayXXf  m(2,2);
  
  // assign some values coefficient by coefficient
  m(0,0) = 1.0; m(0,1) = 2.0;
  m(1,0) = 3.0; m(1,1) = m(0,1) + m(1,0);
  
  // print values to standard output
  cout << m << endl << endl;
 
  // using the comma-initializer is also allowed
  m << 1.0,2.0,
       3.0,4.0;
     
  // print values to standard output
  cout << m << endl;
}

 

3) Addition and subtraction

array + array, array - scalar

int main()
{
  ArrayXXf a(3,3);
  ArrayXXf b(3,3);
  a << 1,2,3,
       4,5,6,
       7,8,9;
  b << 1,2,3,
       1,2,3,
       1,2,3;
       
  // Adding two arrays
  cout << "a + b = " << endl << a + b << endl << endl;
  // Subtracting a scalar from an array
  cout << "a - 2 = " << endl << a - 2 << endl;
}

 

4) Array multiplicaiton

區分於矩陣乘法, 對應位置參數相乘.

int main()
{
  ArrayXXf a(2,2);
  ArrayXXf b(2,2);
  a << 1,2,
       3,4;
  b << 5,6,
       7,8;
  cout << "a * b = " << endl << a * b << endl;
}

 

5) Other coefficient-wise operations

abs(), sqrt(), min() 兩個 array 對應位置取最小

int main()
{
  ArrayXf a = ArrayXf::Random(5);
  a *= 2;
  cout << "a =" << endl 
       << a << endl;
  cout << "a.abs() =" << endl 
       << a.abs() << endl;
  cout << "a.abs().sqrt() =" << endl 
       << a.abs().sqrt() << endl;
  cout << "a.min(a.abs().sqrt()) =" << endl 
       << a.min(a.abs().sqrt()) << endl;
}

 

6) Converting between array and matrix expressions

matrix 有一個 array() 方法, array 有一個 matrix 方法, .array() 和 .matrix() 的返回值可以用作左值或右值.

matrix 有一個 cwiseProduct() 方法, 用於兩個 matrix 對應參數相乘.

Eigen 允許將 array expression 賦值給 matrix 變量.

int main()
{
  MatrixXf m(2,2);
  MatrixXf n(2,2);
  MatrixXf result(2,2);
  m << 1,2,
       3,4;
  n << 5,6,
       7,8;
  result = m * n;
  cout << "-- Matrix m*n: --" << endl << result << endl << endl;
  result = m.array() * n.array();
  cout << "-- Array m*n: --" << endl << result << endl << endl;
  result = m.cwiseProduct(n);
  cout << "-- With cwiseProduct: --" << endl << result << endl << endl;
  result = m.array() + 4;
  cout << "-- Array m + 4: --" << endl << result << endl << endl;
}

稍微復雜點的栗子

int main()
{
  MatrixXf m(2,2);
  MatrixXf n(2,2);
  MatrixXf result(2,2);
  m << 1,2,
       3,4;
  n << 5,6,
       7,8;
  
  result = (m.array() + 4).matrix() * m;
  cout << "-- Combination 1: --" << endl << result << endl << endl;
  result = (m.array() * n.array()).matrix() * m;
  cout << "-- Combination 2: --" << endl << result << endl << endl;
}