Eigen常用操作

        Eigen是一个高层次的C ++库，支持线性代数，矩阵和矢量运算，数值分析及其相关的算法。本例程包含头文件：

#include <iostream>

using namespace std;

#include <ctime>
// Eigen 核心部分
#include <Eigen/Core>
#include <Eigen/Dense>
using namespace Eigen;

#define MATRIX_SIZE 5

1、变量声明

        Eigen 中所有向量和矩阵都是Eigen::Matrix，它是一个模板类。它的前三个参数为：数据类型，行，列。

        1）声明一个2*3的float矩阵

Matrix<float, 2, 3> matrix_23;

        2）内置类型

        Eigen 通过 typedef 提供了许多内置类型，不过底层仍是Eigen::Matrix。

Vector3d v_3d; // 等价于Eigen::Matrix<double, 3, 1>，即三维向量
Matrix3d matrix_33; // 等价于Eigen::Matrix<double, 3, 3>，即三乘三矩阵

        如果不确定矩阵大小，可以使用动态大小的矩阵：

Matrix<double, Dynamic, Dynamic> matrix_dynamic; 
MatrixXd matrix_x; // 更简单点

2、初始化

        Eigen 中所有向量和矩阵都是Eigen::Matrix，它是一个模板类。它的前三个参数为：数据类型，行，列。

Matrix3d::Zero(); // doubel类型，3*3全零矩阵 
Matrix3d::Random(); // doubel类型，3*3随机数矩阵

矩阵赋初值：

Matrix3d matrix_33 = Matrix3d::Zero(); //初始化为3*3全零矩阵

赋值，显示：

// 输入数据（初始化） 
Matrix<float, 2, 3> matrix_23; 
matrix_23 << 1, 2, 3, 4, 5, 6; 

// 输出 
cout << "matrix 2x3 from 1 to 6: \n" << matrix_23 << endl; 

// 用()访问矩阵中的元素 
cout << "print matrix 2x3: " << endl; 
for (int i = 0; i < 2; i++) 
{ 
    for (int j = 0; j < 3; j++) cout << matrix_23(i, j) << "\t"; 
    cout << endl; 
}

        结果：

3、矩阵运算

定义变量：

Matrix<float, 2, 3> matrix_23; // 2*3矩阵，float
matrix_23 << 1, 2, 3, 4, 5, 6;

Vector3d v_3d; // 3*1向量，double
v_3d << 3, 2, 1;

Matrix<float, 3, 1> vd_3d; // 3*1向量，float
vd_3d << 4, 5, 6;

        在Eigen里不能混合两种不同类型的矩阵，这样语法是错误的。

Matrix<double, 2, 1> result_wrong_type = matrix_23 * v_3d;

        同样不能搞错矩阵的维度，这样语法是错误的。

Eigen::Matrix<double, 2, 3> result_wrong_dimension = matrix_23.cast<double>() * v_3d;

        矩阵显示调用：

// 2_3*3_1 = 2_1
Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
cout << "[1,2,3;4,5,6]*[3,2,1]=" << result.transpose() << endl;

Matrix<float, 2, 1> result2 = matrix_23 * vd_3d;
cout << "[1,2,3;4,5,6]*[4,5,6]: " << result2.transpose() << endl;

4、矩阵运算

        一些矩阵运算，四则运算就不演示了，直接用+-*/即可。

Matrix<float, 3, 3> matrix_33;      // 随机数矩阵
matrix_33 << 1, 5, 2, 2, 4, 1, 3, 6, 2;
	
cout << "random matrix: \n" << matrix_33 << endl;
cout << "transpose: \n" << matrix_33.transpose() << endl;      // 转置
cout << "sum: " << matrix_33.sum() << endl;            // 各元素和
cout << "trace: " << matrix_33.trace() << endl;          // 迹
cout << "times 10: \n" << 10 * matrix_33 << endl;               // 数乘
cout << "inverse: \n" << matrix_33.inverse() << endl;        // 逆
cout << "det: " << matrix_33.determinant() << endl;

结果：

5、特征值，对角化

        对角化，求解特征值。

Matrix<double, 3, 3> matrix_33;
matrix_33 <<1, 5, 2, 2, 4, 1, 3, 6, 2;

// 对角化
SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose() * matrix_33);

// 求解特征值
cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;

结果：

6、解方程

        求解 matrix_NN * x = v_Nd 这个方程，N的大小在前边的宏里定义，这里N=50，它由随机数生成。

Matrix<double, MATRIX_SIZE, MATRIX_SIZE> matrix_NN
 	= MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
matrix_NN = matrix_NN * matrix_NN.transpose();  // 保证半正定
Matrix<double, MATRIX_SIZE, 1> v_Nd = MatrixXd::Random(MATRIX_SIZE, 1);

1）直接求逆自然是最直接的，但是求逆运算量大。

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

2）通常用矩阵分解来求，例如QR分解，速度会快很多。

time_stt = clock();
x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
cout << "time of Qr decomposition is "
	<< 1000 * (clock() - time_stt) / (double)CLOCKS_PER_SEC << "ms" << endl;
cout << "x = " << x.transpose() << endl;

  3）对于正定矩阵，还可以用cholesky分解来解方程。

time_stt = clock();
x = matrix_NN.ldlt().solve(v_Nd);
cout << "time of ldlt decomposition is "
<< 1000 * (clock() - time_stt) / (double)CLOCKS_PER_SEC << "ms" << endl;
cout << "x = " << x.transpose() << endl;

结果：

x = -0.177629 0.179189 7.66884e-05 -0.493822 0.568221