eigen库的矩阵操作(Eigen/Dense类)

// 来自高翔SLAM十四讲
#include<iostream>
#include<ctime>
#include<eigen3/Eigen/Core>
// 稠密矩阵的代数运算(逆、特征值等)
#include<eigen3/Eigen/Dense>

using namespace std;

#define MATRIX_SIZE 5

int main(){
	// eigen以矩阵为基本数据元素。它是一个模板类。它的前三个参数为：数据类型，行，列
	// 声明一个2*3的float矩阵
	Eigen::Matrix<float,2,3> matrix23;
	// 同时，eigen通过typedef提供了许多内置类型，不过底层仍旧是eigen::Matrix
	// 例如：Vector3d实质上是Eigen::Matrix<double,3,1>
	Eigen::Vector3d v_3d;
	// Matrix3d实质上是Eigen::Matrix<double,3,3>
	Eigen::Matrix3d matrix33 = Eigen::Matrix3d::Zero();	// 初始化为0
	// 如果大小不确定的矩阵，可以使用动态大小的矩阵
	Eigen::Matrix<double, Eigen::Dynamic, -1> matrix_dynamic;
	// 更简单的
	Eigen::MatrixXd matrix_x;
	// 还有很多类似的类型，自己看吧

	// 下面是对矩阵的操作
	// 输入数据
	matrix23 << 1,2,3,4,5,6;	// 重载了 <<
	// 输出
	cout<<matrix23<<endl;

	// 用()访问矩阵中的元素
	for(int i=0; i<1; i++){
		for(int j=0; j<2; j++){
			cout<<matrix23(i,j)<<endl;
		}
	}

	v_3d << 3,2,1;
	// 矩阵与向量相乘(实际上仍然是矩阵和矩阵)
	// 不同类型的矩阵不能相乘
	// Eigen::Matrix<double, 2, 1> result = matrix23 * v_3d;

	// 应该显式转换
	Eigen::Matrix<double, 2, 1> result = matrix23.cast<double>() * v_3d;
	cout<< result <<endl;

	// 生成随机的矩阵
	matrix33 = Eigen::Matrix3d::Random();
	cout<< matrix33 <<endl;

	// 转置
	cout<< matrix33.transpose() <<endl;
	// 各个元素的和
	cout<< matrix33.sum() <<endl;
	// 迹
	cout<< matrix33.trace() <<endl;
	// 数乘
	cout<< 10*matrix33 <<endl;
	// 逆
	cout<< matrix33.inverse() <<endl;
	// 行列式
	cout<< matrix33.determinant() <<endl;

	// 特征值
	// 是对称矩阵可以保证对角化成功
	Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigen_solver(matrix33.transpose() * matrix33);
	cout<<" eigen value = "<<eigen_solver.eigenvalues() <<endl;	// 特征值
	cout<<"eigen vector = "<<eigen_solver.eigenvectors() <<endl;	// 特征向量

	// 解方程
	// 我们求解matrix_NN * x = v_Nd 这个方程
	// 直接求逆当然是最简单的，但是直接求逆的计算量特别大
	Eigen::Matrix<double, MATRIX_SIZE, MATRIX_SIZE>matrix_NN;
	matrix_NN = Eigen::MatrixXd::Random(MATRIX_SIZE,MATRIX_SIZE);
	Eigen::Matrix<double, MATRIX_SIZE, 1> v_Nd;
	v_Nd = Eigen::MatrixXd::Random(MATRIX_SIZE, 1);

	// 计时
	clock_t time_stt = clock();
	// 直接求逆
	Eigen::Matrix<double, MATRIX_SIZE, 1> x = matrix_NN.reverse() * v_Nd;
	cout<<" time use in normal inverse is "<< 1000 * (clock() - time_stt)/(double)CLOCKS_PER_SEC << "ms"<<endl;

	// 通常用矩阵分解来求，例如QR分解，速度会快很多,然而这边显示时间却慢很多
	time_stt = clock();
	x  = matrix_NN.colPivHouseholderQr().solve(v_Nd);
	cout<<"time use in Qr composition is "<<1000 * (clock() - time_stt)/(double)CLOCKS_PER_SEC<<"ms"<<endl;

	cout<<endl;
	return 0;
}