Eigen 矩陣庫學習筆記

最近為了在C++中使用矩陣運算，簡單學習了一下Eigen矩陣庫。Eigen比Armadillo相對底層一點，但是只需要添加頭文庫即可使用，不使用額外的編譯和安裝過程。

基本定義

• Matrix3f是3*3矩陣，MatrixXf表示矩陣維數不確定，MatrixXf m(3,4)表示3*4矩陣。
• 'MatrixXd是double型，MatrixXf是float`型。

列優先和行優先

Eigen中存儲Matrix用的是column-major，但是初始化賦值的時候是row-major

Matrix3d m;
m << 1,2,3,4,5,6,7,8,9;
	/*	1 2 3
		4 5 6
		7 8 9*/

m(3)=2，而不是4。

矩陣運算

矩陣相乘：m1*m2
element-wise 相乘：m1.cwiseProduct(m2)
沒有除法，但是有倒數，m1.cwiseInverse()，所以m1除m2即m1.cwiseProduct(m2.cwiseInverse())

元素操作

可以對Matrix中的所有元素使用常用函數或自定義函數進行運算，首先需要將Matrix轉換為Array，例如

m1.array().sqrt()

等價於將m1矩陣中的所有元素xsqrt(x)

注意：m1.array()將矩陣m1轉換成了Array，此時m1不再是矩陣，不能夠和其他矩陣出現在同一個表達式中，除非是=賦值運算符。Eigen允許將一個Array賦值給Matrix，也可以將一個Matrix賦值給Array。可以使用.matrix()將Array轉換回矩陣。以下面代碼為例：

out = out.cwiseProduct(bn_var.array().sqrt().cwiseInverse());

上面代碼會報錯，因為bn_var.array().sqrt().cwiseInverse()現在是Array，而不是Matrix。正確的表達式應該是

out = out.cwiseProduct(bn_var.array().sqrt().cwiseInverse().matrix());

可以使用.unaryExpr自定義函數，
出處：https://stackoverflow.com/questions/33786662/apply-function-to-all-eigen-matrix-element @vsoftco

#include <cmath>
#include <iostream>

#include <Eigen/Core>

double Exp(double x) // the functor we want to apply
{
    return std::exp(x);
}

int main()
{
    Eigen::MatrixXd m(2, 2);
    m << 0, 1, 2, 3;
    std::cout << m << std::endl << "becomes: ";
    std::cout << std::endl << m.unaryExpr(&Exp) << std::endl;
}

從CSV中讀取矩陣

函數出處：https://gist.github.com/infusion/43bd2aa421790d5b4582 @Robert Eisele

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

Eigen::MatrixXd readCSV(std::string file, int rows, int cols) {

  std::ifstream in(file);
  
  std::string line;

  int row = 0;
  int col = 0;

  Eigen::MatrixXd res = Eigen::MatrixXd(rows, cols);

  if (in.is_open()) {

    while (std::getline(in, line)) {

      char *ptr = (char *) line.c_str();
      int len = line.length();

      col = 0;

      char *start = ptr;
      for (int i = 0; i < len; i++) {

        if (ptr[i] == ',') {
          res(row, col++) = atof(start);
          start = ptr + i + 1;
        }
      }
      res(row, col) = atof(start);

      row++;
    }

    in.close();
  }
  return res;
}

replicate

vec.replicate(times) 
mat.replicate(vertical_times, horizontal_times)
mat.colwise().replicate(vertical_times, horizontal_times)
mat.rowwise().replicate(vertical_times, horizontal_times)